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Abstract 

In this dissertation I analyze Hamiltonian control of rf-dimensional quantum sys- 
tems as realized in alkali atomic spins. Alkali atoms provide an ideal platform for 
studies of quantum control due to the extreme precision with which the control fields 
are characterized as well as their isolation from their environment. In many cases, 
studies into the control of atomic spins restrict attention to a 2-dimesional subspace 
in order to consider qubit control. The geometry of quantum 2-level systems is much 
simpler than for any larger dimensional Hilbert space, and so control techniques for 
qubit s often are not applicable to larger systems. In reality, atoms have many inter- 
nal levels. It seems a shame to throw away most of our Hilbert space when it could 
in principle be used for encoding information and performing error correction. This 
work develops some of the tools necessary to control these large atomic spins. 

Quantum control theory has some very generic properties that have previously 
been explored in the literature, notably in the work from the Rabitz group. I provide 
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a review of this literature, showing that while the landscape topology of quantum 
control problems is relatively independent of physical platform, different optimization 
techniques are required to find optimal controls depending on the particular control 
task. To this end I have developed two optimal control algorithms for finding unitary 
maps for the problems of: "state preparation" where we require only that a single 
fiducial state us taken to a particular target state and "unitary construction" where 
the entire map is specified. State mapping turns out to be a simple problem to 
solve and is amenable to a gradient search method. This protocol is not feasible for 
the task of finding full unitary maps, but I show how we can weave state mappings 
together to form full unitary maps. This construction of unitary maps is efficient in 
the dimension of the Hilbert space. 

The particular system I have used for demonstrating these control techniques is 
that of alkali atoms, specifically ^^^Cs. The state preparation algorithm was used to 
create a broad range of target states in the 7-dimension F — 3 hyperfine manifold in 
an experiment using a combination of time-dependent magnetic fields and a static 
tensor light shift. The yields from this experiment were in the range of 0.8 — 0.9. I 
have developed another control system for the full hyperfine manifold in the ground- 
electronic state of ^^^Cs, a 16-dimensional Hilbert space, based on applied radio 
frequency and microwave fields. Numerical studies of the state preparation algorithm 
find good operating points commensurate with modest laboratory requirements. This 
system of microwave and rf control also admits a Hamiltonian structure than can be 
used by my protocol for unitary construction. I demonstrate the performance of this 
algorithm by creating a standard set of qudit gates using physically realistic control 
fields, as well as by implementing a simple form of error correction. 
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Chapter 1 
Introduction 



In recent years, it has become increasingly clear that quantum dynamics allow us to 
perform certain tasks in ways that are fundamentally more powerful than their clas- 
sical counterparts. Some examples are quantum computing, quantum cryptography, 
quantum-limited measurement, and many others. This presents a challenge, in that 
quantum systems also appear to be fundamentally more difficult to control. This 
is in part due to the technological challenges of manipulating systems deep within 
the quantum regime, but even with extremely "clean" systems we still have to worry 
about our control routines inadvertently destroying the coherences we are trying to 
protect. Unlike in classical control systems, it is not possible to monitor a quantum 
system passively. This has lead to the study of quantum control theory, the goal 
of which is to develop techniques for implementing non-trivial maps on quantum 
systems in spite of the fragility of quantum states. 

In the work in this dissertation I will be considering what is essentially the "eas- 
iest" classes of quantum control problems. For all of this work, measurement will 
only be considered as a verifier of the control protocol, and as a source for feedback 
control. In addition, for the most part all the states considered will be pure and 
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all the dynamics will be Hamiltonian. Without feedback and open-quantum-system 
dynamics that lead to mixed states the math required to model our control systems 
will be much simpler. I will also only consider two control tasks: state preparation 
and unitary construction. In the first, we would like to find dynamics that arise from 
a physical Hamiltonian that map some particular initial state to an arbitrary but 
fixed final state, and in unitary construction we would like the same dynamics to 
describe a unitary map in its entirety. Even in these idealized conditions, we will find 
that designing optimal quantum control protocols for real physical Hamiltonians is 
a rich and subtle problem. 

1.1 Quantum control 

Quantum control theory comes in two main fiavors: "open-loop" and "closed loop." 
Generically, we have a Hamiltonian for control system of the form 



where we would like to choose the "control waveforms," Cj(t), to implement some 
control task. In open-loop control we must choose the control waveforms without the 
benefit of measurement and feedback. Quantum open-loop control has its origins in 
the fields of physical chemistry and nuclear magnetic resonance (NMR) spectroscopy. 
In physical chemistry the goal is to use laser interactions to drive chemical reactions 
or to excite molecular vibrations and rotations [1, 2, 3]. For NMR imaging, pulses of 
rf magnetic fields are used to produce spin rotations, and by shaping pulses rotations 
can be enacted in a more optimal way [4, 5, 6]. 

Attempts to build a scalable quantum computer have demonstrated the need for 
more accurate quantum control. One of the famous DiVincenzo criteria for imple- 
mentations of quantum computing is the ability to apply elements from a universal 
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set of quantum gates [ ]. While studies in error correction and fault tolerance have 
shown that there is an error threshold below which arbitrary length quantum com- 
puting is possible [ ], the precision of control required to reach this threshold is 
daunting. This has led to many applications of open-loop quantum control in order 
to combat loss of coherence in open quantum systems and to engineer robustness to 
errors in the applied controls [6, 9, 10, 11]. One prime example is dynamical decou- 
pling where one engineers sequences of pulses that prevent loss of coherence between 
a qubit and a non-mar kovian environment to enable quantum memories [12, 13] or 
more recently for protected quantum gates [i^:]. 

Of more relevance to the work in this manuscript is the study of optimal quantum 
control. "Optimal" is a bit of a loaded term in the quantum control theory literature 
and should probably be interpreted according to the colloquial English definition of 
the word. Quantum controls can be optimal with respect to a variety of measures, 
some common examples being the time of a pulse length, fidelity of the time-evolved 
state with the target, or purity of the end product of the control sequence. Any, or all, 
of these measures can be enforced by some sort of cost or objective functional J[cj{t)] 
which must be optimized by a set of control waveforms Cj{t) in order for them to 
be considered optimal. There exist two primary techniques to solve optimal control 
problems. One technique is to solve the problem analytically using geometrical or 
Lie algebraic methods [4, 15, 16]. Alternatively, one can attempt to numerically solve 
these optimization problems either by using gradient search methods [17, 18] or by 
learning algorithms, such as genetic algorithms, where the objective function is cal- 
culated by simply performing an experiment [ ] . I will discuss these two methods in 
more detail in Ch. 2.3 and provide a comparison between the benefits and detriments 
of each method. Framing the problem of quantum control in such a general language, 
optimizing a functional, allows for broad applicabiltiy. Quantum control protocols 
often have elements that are system independent and so the design of new protocols 
for quantum control can impact a wide spectrum applications. Optimal quantum 
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control techniques have been explored on a wide variety of platforms ranging from 
optical [ ], to semiconductors [ ], and superconductors [21, 22]. 

1.2 Atomic spins 

The main application of the control techniques of this dissertation will be towards 
the control of atomic spins. Atomic spins are natural carriers of quantum coherence 
for use in various quantum information processing applications. These systems have 
been of particular interest given their excellent isolation from the environment and 
the available techniques in the "quantum optics toolbox" . Examples include ensem- 
bles of atomic spins as quantum information processing elements [23, 24, 25, 26, 27], 
ion-trap quantum computers [28, 29, 30], and neutral-atom optical lattices [ ]. The 
latter has attracted tremendous attention in recent years, as controllable spin lattices 
are seen as a platform in which to perform quantum simulations of condensed matter 
systems [o^] and studies of topological quantum field theory [ ] . 

Quantum optics is a mature technology. Starting from the work of Glauber, 
Cohen- Tanoudji and others [34, 35, 36], our understanding of the interactions of 
atoms with lasers and other electromagnetic fields has reached an unprecedented 
level. Of particular relevance are the developments in laser cooling and trapping 
[37]. In particular, due to the cesium frequency standard, the atomic properties of 
^^^Cs are extremely well-characterized. The ability to model our control system to 
high precision enables us to start considering quantum optimal control. In the same 
way that liquid state NMR has provided an excellent platform for exploring quantum 
control protocols [4, 6, 5], atomic spins in cold atomic ensembles provide a test-bed 
with unique physical properties that allows for new investigations into control and 
measurement techniques. 

In most theoretical discussions on quantum information theory, the fundamental 
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systems considered are qubits, 2-level quantum systems. This certainly makes sense 
from a theoretical computer science perspective where one can make transparent 
analogies between qubits and classical bits, as well as from an engineering perspec- 
tive since all possible coherent manipulations of 2D quantum systems are simply 
geometric rotations and we would like to easily control our carriers of information. 
From a physics perspective, however, it is not clear whether we should restrict our 
attention to qubits. Atoms have large spins with a rich internal structure. Instead 
of qubits, what if we consider rf-dimensional quantum systems or qudits? This is 
significantly more complicated control problem, however, it allows for the possibility 
of qudits as the fundamental information carriers [ ] as well as the embedding of 
logical qubits in a qudit, which may be advantageous for control or protection from 
errors [ ] . Additionally, manipulating a nontrivial Hilbert space allows us to explore 
interesting dynamics such as quantum chaos [ ^9, 40]. The ability to fully control the 
Hilbert space within the atoms for various applications is an important addition to 
our toolbox of atomic controls. 



1.3 Outline of document 



The theoretical work in this dissertation has, for the most part, been conducted in 
collaboration with the experimental group of Poul Jessen at the College of Optical 
Science, University of Arizona. Many of the results I discuss here have been pub- 
lished previously in three papers. The first. Quantum Control of the Hyperfine Spin 
of a Cs Atom Ensemble [ ], is a direct collaboration with Poul Jessen's lab, in a 
project headed by Souma Chaudhury. For this work I developed a protocol for state 
preparation in atomic spin systems and provided theoretical support for an experi- 
mental implementation of said protocol. The control system in this experiment was 
the magnetic field and AC-Stark shift system initially explored by Silberfarb and 
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Smith [42, 43, 44] and discussed in detail in this thesis. Quantum Control of the 
Hyperfine- Coupled Electron and Nuclear Spins in Alkali Atoms [ ] is a theoretical 
study in which my collaborators and I proposed a new atomic control system that 
uses microwave and rf magnetic fields as the controls. This system should be more 
favorable to implement in the lab and numerical simulations suggest that we can per- 
form state preparation on a space that is twice as big in a time that is about an order 
of magnitude shorter than in the previous control system. Finally, in the last pa- 
per in this dissertation. Constructing Ceneral Unitary Maps from State Preparations 
[40], we developed a protocol to implement the task of unitary construction based on 
our knowledge of how to create good state preparation routines. This construction 
is efficient in the dimension of the Hilbert space of interest and as an example we 
have used this technique to create unitary maps in the microwave and rf magnetic 
field control system. 

I have also participated in several other projects that will not appear in this 
dissertation. Of direct relevance to the contents of this manuscript, there are two 
projects that are nearing completion. The first is a project collaboration with Brian 
Mischuck on the topic of robust control in the microwave and rf magnetic field system. 
The control fields in this system have a geometry similar to the controls in liquid- 
state NMR systems. We are currently working to directly port some of the robust 
control techniques in NMR to this cold atomic spin system. Another project is in 
collaboration with Carlos Riofrio and Steve Flammia regarding state estimation. In 
that project we are trying to understand the power of random unitary dynamics 
with regard to the information content of measurement outcomes states undergoing 
such evolution. In the case where the dynamics describe some fixed orbit in 5u(rf), 
we have proven the system is not driven through an informationally complete set 
of observables. Even though this means there will be density matrices we cannot 
reconstruct perfectly, it appears that on average we can still use this measurement 
procedure to obtain extremely high fidelity estimates for typical quantum states. 
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Table 1.1: Table of published work with location in text. 



The last project I'll mention here is somewhat farther afield and outside the Deutsch 
group. In published work with Dan Browne, Matt Elliot, Steve Flammia, Akimasa 
Miyake and Anthony Short [ ], we were able to show a phase transition in the 
computational power of the cluster state model of quantum computation. This model 
requires a certain quantum state, a cluster state, as a resource for computation. We 
demonstrated that with faulty resource states there is a sharp phase transition in the 
computational power, with respect to the error rate, that occurs at the percolation 
threshold. 

The remainder of this dissertation is as follows. In Ch. 2 I present a background 
review of some of the basics of open-loop quantum control theory. The mains goals 
of this chapter are to understand: when Hamiltonian dynamics are controllable, the 
relative difficulty between unitary construction and state preparation (as explored 
by the Rabitz group [47, 48, 49, 50, 51, 52]), and some practical methods for finding 
optimal control waveforms. In Ch. 3 I describe the physics behind the two atomic 
control systems in this paper. The first control system consists of an "always on" 
nonlinear interaction derived from the AC-Stark effect combined with controllable 
quasi-static magnetic fields. The second system has no laser interaction and instead 
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utilizes magnetic fields oscillating at rf and microwave frequencies. I explain the 
Hamiltonians dynamics of these two systems and rewrite the Hamiltonians in a form 
that is conducive to our quantum control techniques. I will also show under what 
circumstance these systems are controllable. Chapter 4 describes the state prepara- 
tion algorithm I helped to develop. I describe the basic form of the algorithm and 
its application to both control systems, experimentally in the case of the magnetic 
field and Stark shift systems and in numerical simulation for the microwave and rf 
system. Finally, in Ch. 5 I discuss the unitary construction protocol we proposed in 
[^u] and show some examples of unitary matrix construction in the microwave and 
rf control system. 
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Controlling Quantum Systems 



In this dissertation I will be looking exclusively at open-loop techniques for control- 
ling quantum systems. Open-loop control involves designing time-dependent fields 
to generate a dynamical map without using measurement and feedback. This is nice 
in that we are not required to estimate system parameters in real time, but we in- 
stead can perform a more thorough modeling of our quantum system offline. In the 
attempt to find feedback routines one is often forced to consider measurements that 
form a "classical" commutative subalgebra in order to make the problem tractable. 
This is not the case with open-loop control. In some sense, open-loop control allows 
us to explore more of the "quantum'" nature of our control protocols. On the other 
hand, it can be argued that it is really the process of measurement, and in particular 
measurement backaction, that distances quantum control theory from classical con- 
trol theory. With open-loop control, we are really deriving classical control schemes, 
but control schemes for systems that live in complex manifolds such as SU{d). This 
allows us to directly port over some of the work in classical control theory regarding 
control over Lie groups. 

In this chapter I will provide an overview of some very general results from open- 
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loop quantum control, in particular as they apply to the problems of state preparation 
and unitary construction. In Sec. 2.1, I will review what it means to be controllable. 
In the abstract, controllability determines whether a Hamiltonian system has the 
degrees of freedom necessary to perform a given control task. Next, in Sec. 2.2, I will 
look at some of the results on the control landscape topology of state preparation and 
unitary construction. Analyzing the topology allows one to make some surprisingly 
general statements regarding the complexity of finding optimal controls for the two 
different problems. Finally, I'll describe some of the methods we use to find control 
fields. I will discuss two diflFerent classes of optimization algorithms and describe 
some illustrative examples. 

2.1 Controllability 

Before actually trying to control a quantum system, it is a worthwhile endeavor to 
determine whether the system is controllable in principle. When we ignore physical 
constraints like bandwidth and decoherence, what types of control are possible at 
a later, finite time? There are many different aspects of a quantum evolution that 
we might wish to control, and accordingly there are many different concepts of con- 
trollability in the quantum control theory literature. Some common controllability 
questions are whether the available dynamics allow: mappings between arbitrary 
states (pure or mixed), the construction of general unitary maps, or the simulation 
of arbitrary observables. 

In this dissertation, I will primarily consider "unitary controllability", that is 
whether our dynamics allow us to construct any unitary map in a finite time. The 
reasons for considering this type of controllability are twofold. First, unitary con- 
trollability is sufficient for the types of tasks we will consider, i.e. state preparation 
and unitary construction, and indeed most of the pure-state control tasks in the 
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literature. Secondly, the conditions for controllability in the case of unitary control 
are by far the most intuitive and geometrical. 

The conditions for unitary controllability have been studied in depth in the control 
theory literature, [53, 54], and more recently from a quantum information perspec- 
tive [^^]. Formally, we consider a quantum evolution that is governed by a general 
Hamiltonian evolution of the form 



Here, the functions Cj{t) are the control waveforms we are allowed to manipulate. We 
take the operators {i/o, Hi . . . H^} as traceless and Hermitian, which leads to unitary 
dynamics from the group SU (d). In principle, we could consider Hamiltonians with a 
nonzero trace leading to dynamics from but for quantum system global phases 

are irrelevant. For a Hamiltonian system of this form to be considered controllable we 
must show that, starting from the identity operator, we can generate any arbitrary 
unitary operator. More formally, if we let U{t) be the solution of the Schrodinger 
equation 



with U{0) = I, then for some T < oc there exist control waveforms Cj{t) such that 
U{T) = U for any U G SU{d). 

Requiring controllability places constraints on the structure of the independent 
terms in the Hamiltonian {H^^Hi . . . Hn}. In order to generate any element of the Lie 
group SU (rf), it is both necessary and sufficient that the operators {i/o, Hi . . . H^} be 
a generating set for the corresponding Lie algebra 5u(rf). The Lie algebra generated 
by {Hq^Hi. . . Hn} is defined as the closure of the generating set with respect to 
general linear combinations and commutators. 

A Lie algebra is a linear vector space with an algebraic product defined by the 
commutator. We can see that we can generate any linear combination of our initial 




(2.1) 



.dUit) 



H{t)U{t) 



(2.2) 
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set of generators by looking at very short square-pulses according to the Trotter 
formula, where 

Such short pulses are allowed since we assume access to arbitrary control waveforms. 
In addition to linear combinations it is also possible to generate the commutators by 
the approximation 

^-iHiA^-iH2A^iHiA^iH2A _ ^-[Hi,H2]A^ ^ ^2.4) 

The ability to generate, in principle, any linear combination and any commutator 
means that one can simulate any element of the the Lie algebra generated by our 
initial independent Hamiltonians, {i^o, ^^i, • • • , Hn}. 

It is reasonably intuitive to see why {i/g, i^i, . . . , H^} generating 5u{d) will be 
necessary and sufficient for controllability. We can treat the Lie group SU{d) as 
a smooth manifold and Qu{d) as its tangent space. Since we are ignoring physical 
limitations on the control fields Cj{t) we can create infinitesimal displacements along 
the directions described by {i/o, Hi. . . Hn}. To be controllable it is necessary that 
using a finite sequence of these displacements we can simulate a infinitesimal dis- 
placement along any arbitrary direction, since all infinitesimal displacements of the 
identity operator are elements of SU{d). Therefore, it is necessary that the opera- 
tors {i^o, Hi . . . Hn} generate the Lie algebra su{d) through linear combinations and 
commutators. Sufficiency is a consequence of the fact that SU{d) is compact and 
simply connected. This implies that any two elements of SU {d) are linked by a finite 
length geodesic. Access to infinitesimal displacements along all directions in m{d) 
allows us to create an arbitrary geodesic though the identity operator, and thus any 
element of SU{d). 

There are a number of ways to determine whether the independent terms in 
a Hamiltonian control system generate the Lie algebra 5u(rf). The most general 
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approach is to compute the iterated commutators numerically, see appendix B for 
Mathmatica code. We take our initial set of operators and form an orthonormal 



we compute the commutators of all pairs of these Hermitian basis operators and see 
if this results in any operators that have support outside of the initial set. If so, we 
append these to our basis for the algebra. We can look at the commutators of these 
new terms with our basis and iterate until either we span the entirety of 5u(rf), or 
we close on a sub-algebra. 

While this technique, in principle, can work for any set of control Hamiltonians, in 
practice, numerical errors start to become a problem for larger systems. It can also, 
in the worst cases, require calculating something on the order of commutators. 
For large systems it is easier if one can prove controllability analytically by exploiting 
the geometry of the Hamiltonians. 

To close this section I'll prove a simple theorem that we have been able to exploit 
to show controllability in many of the atomic systems I'll be considering in this 
dissertation. 

Theorem 1 In an d- dimensional Hilhert space with d > 2, if one has access to the 
irreducible generators of rotations^ Jx and Jy, then in order to fully control the space 
it is sufficient to add an operator h that has a non-zero overlap (according to the 
trace inner product) with at least one rank- 2 irreducible spherical tensor. That is 



Here we have introduced the orthonormal basis of irreducible spherical tensor 
operators. 



basis with respect to the standard trace inner product 




3 q s± Tr {hT^^^) 7^ ^ { J„ Jy, h}L.A. = 5u(rf). 




m 



(2.5) 



m 
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satisfying the fundamental commutation rules, 

[J„rf] = gTf (2.6) 

where J± = Jx ± iJy. It follows from these commutators that given the set {Jx^Jy^ 
Tg^^} one can simulate any rank- A; irreducible tensor, and since these are an operator 
basis, the generators of rotation can map any rank-/c operator to any other rank-Zc 
operator. With this property we are now prepared to prove a lemma. 



Lemma 1 {/c, J^,Tq } generates Qu{d). 



We prove this by first noting that 

= c,,,Tf + (2.7) 

The exact form of the constants is irrelevant except for the fact that there is always 
some rank- A; tensor for which Ck^q is nonzero. Given this, the proof follows by induc- 
tion. Suppose our library of simulatable operators contains all operators of ranks k 
and k — 1. By commuting some rank k operator with Tq we obtain an operator 
with support on operators of rank k — 1 and A; + 1, thus containing a component in 
the space of rank + 1 operators that is linearly independent from the current set 
of Hamiltonians in our library. Commutation with the generators of rotation allow 
us to simulate all other rank + 1 operators. Since we can simulate all rank-1 from 
the generators { J^^, Jy}^ and the rank-0 operator is the trivial identity operator, it 
follows by induction that we can simulate all rank-/c operators that are supported on 
the Hilbert space, k < d —1. Therefore {Jx-, Jy^T^^} generates 5u(rf).QED 

With this lemma, we see that in order to show theorem 1, we need merely to show 
that the set {Jx^ Jy^h} can simulate the operator Tq^\ We will do this in essentially 



^0 5 
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three steps. Before we start we expand the Hamiltonian h in our spherical basis, 

— 2^k=l 2^q=-k 'k -l-q ' 

d-1 k 

Step 1 : Simulate h = T^^^ + ^Y1 /if ^Tf ) 

k=S q=-k 

To simulate hi we note that h is defined to have some nonzero rank-2 component. 

(2) 

With rotations we can transform the rank-2 component to Tq . Additionally, since 
we have all the rank-1 tensors in our library already, we can remove the rank-1 piece 
of h through linear combinations to yield hi. 

d-1 

Step 2 : Simulate h2 = Tq^^^ + h^^^^T^^^ 

k=3 

Consider the double commutator 

[JM,h]] = '^J2 ^'K^'V- (2-8) 

k=S q=—k 

If we take a linear combination hi — a [Jz^ [Jz-, hi\[ the resulting operator has the same 
coefficients for g = 0. For ^ 0, choosing a = we can sequentially remove all 
rank-2 tensor components, and we are left with /i2. 



Step 3 : Simulate 
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Consider the double commutator 



[Jx^ [Jx) ^2]] 



+ 



k=3 



V(fc-l)A:(A; + l)(fc + 2)(rf) + t['^) 



(2.9) 



We repeat the process in Step 2 to remove the components from /12 with q ^ 



If we now take the hnear combination /12 — 2/12/(^0(^0 + 1)) we remove the T^''°> 

(2) 

component, but are left with a nonzero Tq term. Repeating this procedure for 
/co = 3 . . . (rf — 1) yields an operator that is proportional to Tq^^ This completes our 
proof of theorem 1. 

2.2 Control landscape topology 

In the last section we discussed how to determine whether a Hamiltonian system was 
controllable in principle, but for practical applications we need some way of finding 
the appropriate controls. One would suspect that the relative difficulty of finding 
controls must be very system specific, however, it turns out that it is possible to make 
extremely general statements about the complexity of finding control waveforms. 
This type of analysis derives from studies of the topology of the "quantum control 
landscape" . 



to obtain 




(2.10) 
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Finding optimal quantum controls always corresponds to maximizing some objec- 
tive function J[c] with respect to some control parameters c. Traditionally, J takes 
the form of a fidelity or distance measure and c describes the control waveforms we 
use to drive the system. The quantum control landscape is the multidimensional 
surface described by the value of the objective function as a function of the control 
parameters c. Of particular interest are the critical points on this surface where the 
gradient VcJ = since some must describe the highest quality controls. 

The contents of this section follows from a sequence of papers from the Rabitz 
group on control landscape topology [47, 48, 49, 50, 51, 52]. While the Rabitz group 
has studied a wide variety of control problems, the outcome appears to be the same 
— the landscape topology depends on the dimension of the quantum system and the 
type of objective function, i.e., state preparation, unitary construction, etc., but not 
on any properties of the target or initial states or the particulars of the Hamiltonian, 
excepting controllability. This is an incredibly powerful property since the landscape 
topology alone appears to set the complexity of finding good controls. In this section 
I will paraphrase the arguments in the Rabitz papers in the language I have been 
using in this dissertation. 

The punch-line of this section will be that the problems of state preparation and 
unitary construction have very different control landscape topologies. We will find 
that state preparation, mapping a single initial state to a single target state, has 
an extremely favorable topology that will allow for the construction of very efficient 
search routines for finding control fields. In contrast, the landscape topology of 
unitary construction, mapping the identity to a target unitary operator, is much 
more complex. Numerical surveys [52] suggest that it takes exponentially more 
resources in the dimension of the Hilbert space to search for controls that generate 
unitary maps when compared to those required for state preparation. 

As an aside, in this chapter I will mostly consider the dynamics to be elements 
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of as opposed to SU{d) like in the rest of this manuscript. This assumption 

wiU greatly simplify some of the arguments of this section. Adding a global phase 
is irrelevant to the physics of the problem and in no way diminishes the intuition 
gleaned from these studies. 

2.2.1 Landscape topology of state preparation 

In the problem of state preparation, we would like to map an initially known pure 
state of a rf-dimensional quantum system, to a fixed but arbitrary target pure 
state I/), up to a global phase. The system evolves according to Hamiltonian of 
the form given in Eq. 2.1, and we control this system by specifying the functions 
Cj(t), which are defined from t = [0,T]. For simplicity, instead of using continuous 
functions as our optimization variables, we will assume that the information content 
of the control waveforms Cj{t) can be completely described by some finite length 
control vector c, e.g. square pulses control waveforms or waveforms from cubic splines. 
We can write the Hamiltonian as a function of this control vector, H[c]. From the 
Schrodinger equation, this Hamiltonian leads to a unitary propagator we can write 
as 



where T is the time-ordering operator. Finding good controls amounts to optimizing 
the fidelity between the time-evolved quantum state and the target state, given by 
the objective function. 



The first step to understanding the topology of the quantum control landscape is 
to determine the set of critical points where VcJ[c] = 0. A perfect state preparation, 
J = 1, is an extremal point of the control landscape and thus must be a member of 



U[c] =r [e-^^o^^'^i^i 



(2.11) 




(2.12) 
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the set of critical points. To simplify our calculation of the critical points we define 
the Hermitian matrix A as the logarithm of [/, U[c\ = e~^^[^l. This will always exist 
since U is unitary, however, A is in general an extremely complicated functional of 
c. For the following calculation I'll write the eigen-decomposition of A as 

A = J2('Mj){<^j\- (2-13) 

J 

An alternative description of A is in terms of some orthonormal Hermitian basis, E^^ 
so that 

A = J2AvE,, (2-14) 

T] 

with Aj^ = Tt{AEj^). The canonical basis we will use consists of (f terms of the form 
+ |e/c)(ej|) /\/2 for j < A; and {-i\ej){ek\ + i\ek){ej\) /V2, also for 

j < k. 

The key insight from ['^^] is that it is possible to remove essentially all of the 
particulars of the Hamiltonian dynamics from the condition = by a very 

simple argument from controllability. We can use the chain rule to rewrite Vc-^^ = 
as 

k k rj ' rj ' 

It follows from controllability that the vectors VcAj^ are linearly independent for 
different ry, and so the derivatives dJ /dAj^ must each independently go to zero. This 
leaves us with new constraint equations that take the form 



which has removed all the dependence on H and c. 



To see that the vectors Vc^ry are linearly independent, consider the following ar- 
gument. One of the necessary implications of controllability is that we can construct 
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any unitary matrix of the form exp(— iaE"^) for all /i and a = [0, oc). This means 
there must always be a control vector c such that Aj^[c] = a and A^^^[c] = 0. We 
can now show why the gradient vectors must be independent through a linearity 
argument. Assume that there exist coefficients such that 

Ve^^ = J]^,VeA,. (2.17) 

This implies that 

= Vc(A^-5]M,) (2.18) 

or 

A^^Y.P,A, + C, (2.19) 

where C is a constant with respect to c. In this case the only unitary matrix we 
can construct of the form exp(— zaE^^) is a single matrix, exp(— z(7£^^). This implies 
that any system where the vectors Vc^ry are linearly dependent is not controllable, 
and so by the contrapositive, if our system is controllable, the vectors Vc^r^ must be 
linearly independent. 

We can look at dJ/dAy^ in more detail by first evaluating dU/dAy^. We do this 
by explicitly differentiating the operator and expressing A in its eigenbasis 

dU 

OA, 



^^^-^A^l-s)^A^-^As 

dA^ 



Jo 

= / ds\am)e-"'-^'-'Ham\Er,\an)e-"'-'{an\ 

= ^\am){0'm\Erj\an){an\F{a,ri,an)- (2.20) 
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Here F{am^ the result of the integral and has the form 



, (—i)e '""^ if 

Fiam,a^)={ \_1 _ " (2.21) 

it am T 



We can plug this back into Eq. 2.16 to get the set of constraint equations 
n - — 

= J]^(«|f/^|/)(/|am)(a^|£;^|a„)(a„|i)F(a^,a„) + c.c. 



(2.22) 



Since these equations must be zero for all ry, they must also be zero for any general 
linear combination. In particular, by making a unitary transformation, we obtain 
that for all r and s 



= ^(a,|E^| 



dJ 
"dA, 



^^{i\U''\f){f\am){am\E^\an){as\E^\ar){an\i)F{am,an) + c.c. 

(2.23) 



One of the consequences of our choice of basis is that it is easy to show that 

J](^l|E,|0l)((/)2|^,|^2) = (^l|^2)(02|</>l), (2.24) 

which leaves us with 

= J^(i|t/^|/)(/|a^)(a^|a^)(a5|a^)(a^|i)F(a^,a^) + c.c. 

= (^|t/t|/)(/|a,)(a,|i)F(a„a,) + c.c. (2.25) 
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To simplify this expression further by we write = U\i). When we remove \i) we 
are left with 

= {a,\U^i^){i^\f){f\ar)F{ar,a,) + c.c. 

= (a,|^)(V'|/)(/|a,)e^"^F(a„a,)+c.c. (2.26) 

At this point it helps to look separately at the cases where = and ^ 
in order to see what restrictions are placed on the time-evolved state, |?/;), by these 
equations. When = the constraint equations reduce to 

= (a,|^)(t/.|/)(/|a,)e-'-(-i)e--'- + c.c. 

= -z(a^|^)(^|/)(/|a^) +C.C. 

= -i{{ar\ij){m{fK) - {ar\f){m{i^\ar)) 

= -i{{ar\m{i^l\f){f\]\ar)). (2.27) 

The equations concerning 7^ are a bit more tricky to deal with. We first 
define, only for ^ a^, the function 

^-i(ar-as) _ 1 

G{ar, a,) = e^"^F(a^, a,) = . (2.28) 

G has two properties of note that one can easily show: Re(G') 7^ and G*{ar, a^) = 
—G{as, ttr). We can rewrite the constraint equations for a^. ^ in terms of G as 

= {a,\^) {ij\f) {f\ar)G{ar, tts) + c.c. 

= (a,|V)(^|/)(/|a,)G'(a„a,) + (a,|/)(/|V)(^|a,)G'*(a„ a,) 
= {as\tp){^p\f){f\ar)G{ar,as) - {ar\f){f\tp){'ip\as)G{as,ar). 

(2.29) 
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While this doesn't immediately look to be an improvement we can look at the sum 
of the two constraint equations for (r, s) and {s, r) to obtain 

= {as\xp){i>\f){f\ar)G{ar,as) - {ar\f){f\ij){ij\as)G{as,ar) 

+ {ar\ip){ij\f){f\as)G{as,ar) - {as\f){f\ij){ip\ar)G{ar,as) 

= (a,|[|^)(^|,|/)(/|]K)G'(a,,a,) + (a,|[|^)(^|,|/)(/|]K)G(a„a,) 

= (a,|[|tA)(V'|,|/)(/|]|ar)G(a.,a,) + (a,|[|V')(V^|,|/)(/|]|a,)G*(a„a,) 

= {asm{^P\,\f){f\]\ar)iG{ar,as) + G*iar,a,)). (2.30) 

Remembering that the real part of G is never zero, we can combine this result with 
the outcome of Eq.2.27 to determine that if VcJ = 0, 

[l^)(^U/)(/|]=0. (2.31) 

That is, the time-evolved state must commute with the target state. The implication 
is that either, the time-evolved state is orthogonal to the target state, or, up to a 
global phase, it is equivalent to the target state. Therefore, when VcJ = 0, J = 0, 1. 
This result is independent of the target and initial states as well as the details of the 
Hamiltonian evolution. 

This implication that VcJ = if and only if J = 0, 1 dramatically impacts the 
ease of search when looking for optimal controls. The control landscape has no sub- 
optimal traps. It isn't necessary to resort to complicated algorithms like genetic 
searches or annealing methods to find global optima. Instead, local algorithms, like 
gradient searches, should converge on globally optimal controls. 

By analyzing the topology of the set of critical points J = 1, we find the structure 
of state preparation is even more favorable. It turns out the set of good controls form 
a manifold. We can see this by looking at the set of unitary operators W C SU {d) 
for which = 1. Since SU{d) is invariant to right multiplication we can 
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re-express the critical set as >V = VR, where R G SU{d) and satisfies R\i) = |/). 
Now the condition for J = 1 is = 1, implying that the only requirement 

on V is that its elements have |/) as an eigenstate. The elements of V are allowed 
to have any unitary action on the orthocomplement of |/), and so V is isomorphic 
to U{d — 1). This is U{d — 1) and not SU{d — 1) due to the unconstrained phase 
associated with the eigenvalue of |/). 

The importance of the critical points, J = 1, forming a smooth submanifold of 
U [d) is that it lends the state preparation problem a certain amount of robustness to 
variations in the control fields. The optimal control fields form a large plateau in the 
control landscape as opposed the case where high fidelity controls could have been 
represented by isolated points in the landscape. When we perturb the control fields, 
only the resulting displacements in SU (d) that have support outside of the tangent 
space of the critical submanifold will lead to a decrease in fidelity. The dimension 
oiU{d— 1) is [d — 1)^, which is a very large fraction of d^ — 1, dimension of SU{d). 
The difference between these two dimensions is only 2g? — 2, which should come as 
no surprise since it is the exact number of parameters necessary to describe a pure 
state. 

2.2.2 Landscape topology of unitary construction 

In the problem of unitary construction, instead of solely mapping one known state 
to some other state, we would like the final, time-evolved unitary map, t/[c], to be 
some particular, but arbitrary, unitary map V G SU{d). We can quantify how close 
the time-evoloved unitary is to the target by the Hilbert-Schmidt distance 



\\U[c]-V\\hs 



VTr|t/[c] - y|2 = V2rf - 2Re {Tv{VW[c])) 



(2.32) 



from which we obtain 



J[c] = 2Re {Tx{V^U[c\)) , 



(2.33) 
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as the objective function we would like to maximize for perfect unitary construction. 

The analysis of this problem proceeds very similarly to that of state preparation, 
even though the two objective functions are quite different. We would like to deter- 
mine the nature of the critical manifolds for which VcJ = 0. First, using the same 
decomposition and insights on the nature of controllability as we did in Eq. 2.15, 
we remove all dependence on the particulars of the evolution to get the independent 
constraints 

Tr( ) +C.C. = V77. (2.34) 



We can directly plug in the value of dU[c\/dAjj from Eq. 2.20 into the set of 
constraint equations to obtain 

^{am\V^\an){am\Erj\an)F{am,an) + C.C. = 0. (2.35) 

To get this into a more manageable form we make the same change of basis as in 
Eq. 2.25 yielding 

dJ 

= ^{^s\Erj\ar)-Q^ = {as\V^ar)F{ar, as) + c.c. (2.36) 

r] ^ 

Again, we look separatly at the cases = and 7^ a^, but this time we will 
first look at the case 7^ a^. We can explicitly write out F(a^, a^) and simplify to 
get 

= {as\V^ar) F{ar^ as) + c.c. 

g ICLr g ICLs 

= (a^iy^'la^) h c.c. 

j» ft s 

= ^ ((a,|0f/|a,) - (a,|t/0|a,)) + c.c. 

a J' a s 

= -^—{as\[V\U]\ar)+c.c. (2.37) 

a^ a q 
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Since these equations must be true for all ^ ag we have that both the real and 
imaginary parts of the off-diagonal elements of [V^ ^ U] must be zero. The diagonal 
components of the commutator must be zero independently in this particular basis 
since |a^) is an eigenstate of U and so 

{ar\[V\U]\ar) = {ar\V^U\ar) - {ar\UV^ar) 

= {ar\V^ar)e-'''^ - e-^^^(a^|y Vr) = 0. (2.38) 



These constraints lead to a similar commutator restriction as in Eq.2.31. That 
is, when VcJ 0, 

[V, U] = 0. (2.39) 

Unlike in state preparation this is not the whole story. In the state preparation 
problem both the evolved and target states were rank-1 projectors, and since global 
phases are irrelevant, the map was defined solely by its eigenvectors. In order to 
construct a full unitary map we must not only consider the eigenvectors of the evolved 
operator, but also their eigenvalues. For this we need the equations corresponding 
to Qr = ag. We know that U and V have simultaneous eigenstates, and write the 
eigenvalues of V as (as|]/|as) = e~^^^ This leads to the constraint equations 

= {as\V^as)F{as, as) + c.c. = e^^^(-z)e-^^^ + c.c. = 2 sin (6, - a,). (2.40) 

For this expression to be zero, bg — ag = n^yr, where Ug is an integer. 

We can now consider what VcJ = implies about the value of J. If the gradient 
of J is zero then 

J[c] = 2Re {Tt{V^U)) = 2Re e-^(«^-^^) j = 2 J^(-l)^^ (2.41) 

This leads to d different values of the objective function ranging from — 2rf, —2d + 
4, . . . , 2rf — 4, 2(i, or Hilbert-Schmidt distances 0, 2, . . . 2d. Unlike in the case of state 
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preparation, when optimizing unitary maps there are d + 1 critical manifolds for 
which the gradient is zero. 

In order to more fully understand the topology of the control landscape for this 
problem, we can look at the group structure of the critical manifolds exactly like in 
[51], Since U{d) is invariant under left multiplication we can make a transformation 
to some W G U{d) such that W = V^U. Under this mapping, the subspace where 
ReTr(VI^) is equal one of the critical values of J is topologically equivalent to the 
subspace where ReTT{V^U) is equal to the same critical value. In one of these critical 
manifold, the matrix elements of W have the form 



W has a block structure of the form Id-n © —In, where n is the number of eigenvalues 
with value —1. In fact, the critical manifold is all such W G U{d) that have this 
eigenspectrum since ReTT{TWT^) = ReTr(VI^) for all T G [/ (d). From here on we will 
label the separate critical manifolds by a canonical representative Wn that is diagonal 
and whose matrix values on the diagonal are arranged (1,...,!,— 1,...,— 1). 

More formally, the set Orb(W^) = {TW^Tt : T G U{d)} is defined as the orbit 
of the group action of U{d) with respect to Wn- Since U{d) is a compact Lie group, 
the orbits form smooth submanifolds of U{d). Additionally, while Orb(W^) is not 
necessarily a group, it is diffeomorphic to the quotient group [/(rf)/Stab(W^). Here 
Stab(W^) is the stabilizer group of Wn in defined as Stab(W^) = {R G U{d) : 

RWnR^ = Wn}- Because Wn has the block structure Id-n © —^m Stab(W^) is simply 
Stab(W^) = {Ud-n®Un : Ud-n ^ U{d — n)^ Un G U{n)}. The critical submanifold has 
the structure of the Grassmannian manifold, that is the manifold of U{n) subspaces 
of U (d) or 



(2.42) 



i i 



G{n, d) 



U{d) 



(2.43) 



U{n) X U{d-n)' 
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The dimensionality of these manifolds is 

dim(G(n, d)) = dim(C/(d)) - (dim(C/(d - n)) + dim(C/(n))) 

= (f-{{d-nf + n^)^2n{d-n), (2.44) 

Unlike in state preparation, where the optimal critical submanifold had a rela- 
tively high dimension, for unitary construction the optimum is a single point. In 
this control landscape, it is the suboptimal manifolds that have dimensions on the 
order of rf^. If any of the suboptimal manifolds were traps, using local searches 
would become hopeless. We can examine the curvature in the vicinity of the critical 
manifolds to determine whether they are saddles or local maxima by computing the 
Hessian. This wasn't necessary in the case of state preparation since the only critical 
manifolds were at the extrema, and thus had to be either maxima or minima. 

The Hessian is essentially the second derivative of the control landscape and has 
matrix elements defined by 

= 

The eigenvalues of the Hessian matrix describe the curvature of the control landscape. 
The key quantity of interest is the sign of the eigenvalues, which determine whether 
the suboptimal manifolds are traps or saddles. Like the rest of the analysis of the 
landscape topology we'll look at variations with respect to the manifold of unitary 
operators as opposed to variations in the control fields. 

The easiest way to understand the eigenspectrum of Ti is to look at the Hessian 
quadratic form. We can rewrite our objective function J[c] as a functional of the 
time evolved unitary map, [/, 

J[U] =2ReTr(yt[/). (2.46) 
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The Hessian quadratic form, Hh{U), is the second order term of the Taylor expansion 
about h of J[e~^^[/], where h is an arbitrary infinitesimal Hermitian operator. The 
Taylor expansion up to second order of J is 

J[e-'^U] = 2ReTr {V\l - ih - h^/2)U) . (2.47) 

Therefore, the Hessian quadratic form is 

HhiU) = -ReTr (V^h'^U) . (2.48) 

We would like to evaluate this quantity when Un is a member of one of the critical 
submanifolds. If we write the matrix values of h in the eigenbasis of V as {Pj\h\f]j) — 
"jjj and {Pj\h\(3k) — 7jk + ^Vjk^ where the 7's and ry's are real, we are left with, 

HhiUn) = -ReY,{Ps\UV^h^m 

s 

= -ReJ2{-ir^{/3,\h'\Ps) 

s 

s,t s>t 

(2.49) 

The independent terms in this sum give us the eigenvalues of H. We can enumerate 
the number of positive, 7"^+, negative, H-^ and zero, T^o, terms in this sum to obtain 

n+ = n^, H-^{d-nf, no^2n{d-n). (2.50) 

The size of the zero eigenspaces confirm our previous geometric arguments. 

The eigenspectrum of the Hessian tells us that the topology has no traps, only 
saddles. Ruling out the possibility of traps might give us hope that the same local 
searches that are efficient in the problem of state preparation should apply here. That 
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is not the case. Numerical simulations have shown [52] that the amount of compu- 
tational resources necessary to optimize unitary maps grows exponentially with the 
dimension of the system. It is not fully understood why the resources should scale 
exponentially with this topology. One clue that is suggested from the numerical stud- 
ies of the landscape is that the path traversed by the optimization increases linearly 
with problem size for optimizing full unitary matrices, while with state preparation 
this distance is roughly constant. For the problem of state preparation, any arbi- 
trary control vector is close to some optimal control. This is impossible in the case 
of unitary construction when the optimal control is a solitary point. 

2.3 Generating optimal control waveforms 

We have discussed how to determine whether a Hamiltonian system is controllable 
and the relative difficulty of the two types of control tasks in this dissertation. In 
this section I'll review some of the techniques for the practical construction of control 
waveforms. For the most part, the algorithms used to construct controls fall into 
one of two broad categories which I will label "stochastic searches" and "geometric 
constructions." In this section I will discuss the relative strengths and weaknesses of 
these two approaches and describe some of the representative algorithms from each 
set. 

2.3.1 Stochastic searches 

The algorithms that I will refer to as stochastic search algorithms all involve the same 
basic steps. First, we select an arbitrary control field from some distribution to serve 
as a random seed. We then use this seed to perform an optimization that attempts 
to maximize our objective function. If this optimization yields controls that are in- 
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sufficient for our needs, we simply draw a new random seed and repeat the process. 
Eventually, this process will find control waveforms such that the value of the objec- 
tive function is arbitrarily close to the global optima. Some optimization routines 
such as simulated annealing or genetic algorithms incorporate the stochasticity in a 
more regular way, but the end result is the same. 

This kind of technique represents a brute force approach to finding optimal con- 
trols. We essentially ignore everything we know about the underlying physics of the 
system and make random guesses that we hope are in the neighborhood of a global 
optima or at least a path to a global optima. Ignoring the structure of the prob- 
lem comes at a steep price. For some problems the time required for these types of 
algorithms to converge on an acceptable answer may become prohibitive, e.g. the 
computational complexity scales exponentially in d. 

While it may seem silly to try to guess the answer, the fact that we can ignore 
all of the particulars of a problem is also a virtue. These types of optimization 
procedures can be constructed for any type of control problem. Stochastic searches 
always represent a possible avenue of last resort, and for small dimensional problems 
the asymptotic scaling can be insignificant. Also, from a practical perspective, since 
these algorithms are all very similar, once one has implemented a stochastic search 
algorithm for one problem, it is almost trivial to retool it for use on a different 
physical system. The ease of implementation is furthered by the availability of canned 
numerical solvers for these search problems for most computer algebra packages. 

Stochastic search algorithms become important when we consider the results from 
Ch. 2.2 regarding the landscape topology of state preparation. State preparation 
has a topology that is extremely favorable towards stochastic searches since it has 
no suboptimal traps and the optimal points form a submanifold of reasonably high 
dimension. With the problem of state preparation we can be sure that a random 
guess not only will always lead us to a global optima but also will be able to do 
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so for local searches. To solve a state preparation problem we do not need genetic 
or simulated annealing algorithms, but instead can make do with simpler gradient 
ascent techniques. For this reason gradient searches have yielded some very powerful 
optimal control search routines. 

Gradient searches are most simply explained in a couple lines of pseudocode. 
c = RANDOM 
while ||VcJ[c]|| > 5 
c + eVcJ[c] 

end 

output c 

We start from a random seed and calculate the gradient of J. As long as we are not 
at a critical point already, the algorithm takes a small step in the direction of the 
gradient. If e is small enough the algorithm will converge on a critical point where 
Vc J[c] = 0. In the problem of state preparation this will always be a global optima. 
There are extra bells and whistles one can add to the algorithm, e.g. adaptively 
choosing e or adding some stochasticity to the objective to help traverse saddles, but 
gradient searches will still find global optima reliably for only the most simple topolo- 
gies. Luckily for us, state preparation has such a topology. For unitary construction, 
we must consider different methods for all but the smallest size problems. 

2.3.2 Geometric constructions 

The algorithms for generating quantum controls that I have described as geometric 
constructions are many and varied. Depending on the structure of the Hamiltonian 
and the type of control problem one is considering, it is occasionally possible to find 
deterministic algorithms that create good control waveforms. These constructions 
are particularly nice since they generally require only minimal computational re- 
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sources, e.g. solving a simple geodesic equation [1]. While we know that it is easy 
to construct state preparations using stochastic searches, geometric constructions 
have, until very recently, been the only way to construct unitary operators with a 
reasonable asymptotic scaling. 

The limitations with these approaches is that SU (d) is a pretty complicated place. 
Unlike the broad applicability of stochastic techniques, the set of problems for which 
we understand the geometry well enough to develop efficient unitary constructions 
is limited. Additionally, geometric controls very often aren't optimal with respect 
to measures such as the total time of the control waveform or the robustness to 
errors. When performing a stochastic search we could simply make adjustments to 
the objective function, but with a geometric construction, altering the objective can 
very easily destroy the geometric property one is exploiting. 

Perhaps the simplest type of geometric construction for unitary matrices is that 
of the Euler angle construction for 2-level systems. While a trivial example, but 
it does encapsulate some of the flavor of these techniques. Our understanding on 
how to construct a 2-level unitary matrix relies on the fact that SU (2) is a double 
cover of S0{3)^ the symmetry of the 2-sphere, which is geometry about which we 
understand well. Given two Hamiltonians, Hq and i^i, we can flnd a set 
trivially, using only trigonometric functions, such that U = e~^<^^oe~^f^^^e~^^^^ ^ for 
any U G SU{2). There does, however, most likely exist some continuous control 
waveform Co{t)Ho + Ci{t)Hi that creates this transformation with a smaller energy 
cost. 

Of course, the main limitation of the Euler angle approach is that it fails for any- 
thing other than 2-dimensional systems. The special unitary group is only isomorphic 
to a sphere for = 2. Furthermore, while there exist some similar constructions in 
higher dimensions, e.g. the Cartan decomposition in = 4, these decompositions 
place requirements on the nature of the Hamiltonian beyond simple controllability. It 
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is more interesting to look at families of geometric constructions that are applicable 
to any dimension. Since there is really no overarching algorithm that describes all 
geometric constructions, I will describe two particular examples from the literature 
that exemplify some of the powers and limitations of this approach. 

Unitary construction from a QR decomposition 

A procedure to exactly construct general unitary operators on a qudit was put for- 
ward in [ ]. This construction requires some very specific Hamiltonian structure. 
The Hamiltonians all come in pairs and these provide controUably on a 2d subspace 
of the form 

= \j){k\ + \k){j\, = -t\j){k\+t\k){j\. (2.51) 

We can define a coupling graph for this system as a graph where the vertices are 
the basis states of our qudit and the edges connect the coupled 2d subspaces. It 
is possible to show the system is controllable if and only if this coupling graph is 
connected. 

An arbitrary unitary map on this system can be implemented through a method 
that is derived from the QR decomposition. All invertible matrices V can be written 
in the form V = QR where Q is a sequence of Given's rotations, Q = G1G2 . . . 
and R is upper triangular. If is a unitary matrix, R must additionally be diagonal. 
A Given's rotation is rotation in a plane spanned by two coordinate axes, i.e., 

G = I+{cos9- + .(cos^ - l)\k){k\ + {sm9)\j){k\ + {sm9)\k){j\. (2.52) 

This decomposition provides a method to construct a general unitary matrix by 
way of backwards-evolving the target to the identity. We simply find a sequence 
of rotations in our 2d subspaces such that G^^ . . .G\V is diagonal. We can do this 
by finding rotations that sequentially set the off-diagonal matrix elements of V to 
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Figure 2.1: In (A) we have an example of a system that allows for the unitary 
construction technique in this section [ ]. The system of interest is the electronic 
ground state of ^''Rb. The coupling Hamiltonians are realized by two lasers driving 
Raman transitions. These resonances can only couple hyperfine levels satisfying the 
selection rules A^^ = 0, ±1, which leads to the connected coupling graph in (B). 



zero. There is a systematic way to set these elements to zero using spanning trees 
of the coupling graph. Details can be found in [15]. Once we have a diagonal 

(z) 

matrix it is simple to remove the phases by considering rotations along Hj in our 
2-dimensional subspaces. We can create these easily enough since the 2d subspaces 
are fully controllable. Now that we have a construction for and we can simply 
apply the time-reversed fields to map the identity to V. 

It should be noted that not only does this technique only work for a very restricted 
class of control Hamiltonians such as the one in Fig. 2.1. This construction does not 
make particularly efficient usage of the available resources. One can discard couplings 
terms and as long as the graph remains connected it turns out that total time of the 
construction remains constant. This construction is more of the form of a proof of 
principle, similar to the Trotter expansion from Ch. 2.1 in that the construction is a 
sequence of single-Hamiltonian propagators. Unlike the Trotter expansion this has 
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no infinitesimals and thus could be used in practice. While this construction is not 
time-optimal, the length of the waveforms is still scales polynomially in rf, and more 
importantly describes a deterministic algorithm. 

Time-optimal control with Riemannian symmetric subspaces 

In [4], the authors developed a very clever way to find time-optimal controls for 
certain types of spin systems by solving a simple geodesic equation. The geomet- 
ric requirements for this scheme are that we have a standard control system on 
a rf-dimensional Hilbert space, given in Eq. 2.1, where the constant term is much 
weaker than than the time-dependent pieces, ||i^o|| <^ Ikj^jll- Furthermore, the 
system must contain of Riemannian symmetric subspace which has the following 
form. We will label the Lie algebra generated by just the time-dependent terms, 
{i/i, . . . ^n}, as fi, with corresponding Lie group K. The (right) coset space 
of the respective Lie group, SU{d)/K^ must be a Riemannian symmetric subspace. 
More precisely, let m denote the orthogonal complement of fi in 5u(rf). The coset space 
SU (d) /K 18 Riemannian symmetric if all elements in fi and m satisfy the commutator 
relations 



This is obviously a fairly restrictive property. One common example however is in 
SU{4:) where K = SU{2) x SU{2). That is, the drift term describes a coupling term 
between two qubits and we completely control the single qubit Hamiltonians. 

The importance of this type of system is that there is now an equivalence between 
finding controls that minimize the time T such that 



C t [t,m] c m [m,m] c fi. 



(2.53) 



U{t) 



iH{t)U{t) 



U{0) = I, U{T) = V 



(2.54) 



and finding time-optimal controls X such that 



Pit) = X{t)Pit) 



P(0) = I P{T) = KV. 



(2.55) 
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Start 



Figure 2.2: Figure from ['] that schematically describes the time optimal controls. 
The dotted line denotes the optimal control to drive the system to V in the shortest 
time. Since movement along the cosets has essentially no cost the algorithm mini- 
mizes how long we must evolve according to the drift term, which is the only way to 
move between cosets. 



Here X belongs not to the entire unitary group, but simply X = Adx(i^o) = 
{k~^Hok : k G K}. This second optimization is much easier because since the 
solution basically describes geodesic equation. 

Instead of moving on SU{d) the second optimization moves through cosets, see 
Fig. 2.2. Since the time-dependent terms are much stronger than the drift term, 
moving within a coset has essentially no cost. Our optimizations simply needs to 
find the point on our current coset where Hq describes the greatest rate of change. 
This means we can use a simple greedy search to find time optimal controls since 
we optimize that rate of change independently at each point. It is crucial that the 
coset space is Riemannian symetric since otherwise the optimal controls may involve 
backtracking, which makes a greedy search impossible. 

This construction is very nice in that it provides not only the optimal controls 
with respect to the Hilbert-Schmidt norm, but also the optimal controls with respect 
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to the duration of the control pulses. The final algorithm for constructing controls 
is simple and deterministic. Again, however, the restrictions on the character of the 
control Hamiltonians reduce its applicability to a small set of physical systems. 
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Alkali Atomic Systems 



As stated in the introduction, atomic spins are a natural system to consider for 
storing and manipulating quantum information. Because of the advances in laser 
cooling, ensembles of alkali atoms are a natural system to explore. When the atoms 
are cold, their motion is negligible and they can be considered to be frozen in space 
over the time scale of interaction. This vastly simplifies the description and allows 
us to focus solely on the internal dynamics. The internal state of alkali atoms is 
dependent only on a single valence electron plus nuclear spin, leading to a hydrogen- 
like level structure. For many isotopes this leads to electronic ground states that 
have a non-trivial number of hyperfine states, e.g. ^^^Cs has a nuclear spin / = 7/2 a 
thus 2(2/+ 1) = 16 sublevels. Since these atoms are neutral and have no dipole in the 
ground state, they are extremely well-isolated from the environment. Furthermore, 
we have easy access to the mature technology of diode lasers that can be tuned to 
the Dl and D2 resonance lines in alkalis, which lie in the near infrared. 

We seek to control the quantum state of a multilevel atom. Though single- 
atom addressing and measurement are possible [56, 57, 58], in practice we consider 
ensembles of uncorrelated particles. To the degree that the atoms are identically 
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Figure 3.1: The level structure of ^^^Cs for the Dl and D2 line. Our control system 
of interest is the 5^1/2 ground state, highlighted in orange. 



prepared and uniformly addressed, with no interactions between them either from 
interatomic forces or through measurement backaction, we can take the joint state of 
the system as effectively identical copies, p^^ . More general many-body control is 
not considered here. Restricting then to a single atom, the relevant Hilbert space of 
an alkali atom in its electronic ground state is the tensor product space of electronic 
spin S and nuclear spin / subsystems, H = l)s®^i- Given the single valence electron 
S' = 1/2, the Hilbert space is spanned by two irreducible subspaces of total angular 
momentum F± = / ± 1/2, such that 7"^ = [)+ © [)_. With -^^^Cs, where the nuclear 
spin is 7/2, these spin manifolds are = 4 and F_ = 3, see Fig. 3.1. 

The Hamiltonian describing the atom and its interactions with external magnetic 
and electric fields in the electronic ground state is given by 

H = //atom + Hmag + //laser = AI - S - p • B{t) - ^E,{tyEj{t)a,j. (3.1) 

Throughout this discussion I will set h = 1. For all the work in this dissertation the 
dominant term in this Hamiltonian will be the hyperfine interaction, ^I-S. In units of 
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Plank's constant, in cesium the strength of the hyperfine couphng is A/h = 9.2 GHz. 
The strength of the apphed magnetic fields will at most be 2/x^5o//i ^ 1 MHz for 
a static bias field but will more typically be on the order of 2/x^|B(t)|//i ^ 10 kHz 
for our time-dependent control fields. The goal of this chapter is to rewrite this 
Hamiltonian in a way that is conducive to the types of control techniques we discussed 
in the previous chapter, as well as showing the resultant systems are controllable. 



3.1 Quasi-static magnetic fields and light shift 

One approach to controlling atomic spins is with Zeeman and AC-Stark shift inter- 
actions [42, 43, 44]. In this control system, the space of interest is restricted to the 
manifold F_ . For the remainder of this section I will label the irreducible generators 
of angular momentum on this space as simply F. Restricting to F, we can write 
H = FfHFf. Since the hyperfine interaction 

• S = ^ - /2 - S^) , (3.2) 

is a constant when reduced to one spin manifold we are left with two terms in our 
control Hamiltonian 

= -/X . B(t) - ^-E,{tYE^{t)a,,. (3.3) 

The nuclear magneton /x/ is about three orders of magnitude smaller than the 
Bohr magneton /x^. We can thus, with high accuracy, write the magnetic field 
Hamiltonian as an operator purely on the electronic spin 

^MAG = 2/XBB(t) . S. (3.4) 

In the linear Zeeman regime, with no resonant eflFects, this Hamiltonian approxi- 
mately preserves F and can be written according to the Lande projection theorem 
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HuAG- f^BgfB{t)-F. (3.5) 

In the experiment I will discuss in Chapter 4, we controlled just the x and y compo- 
nents of the magnetic field. We can combine constants to get the Larmor frequencies, 
in the two directions to write 

^MAG = ^At)F. + ^y{t)Fy (3.6) 

It should be clear that magnetic fields only generate rotations, and thus a rep- 
resentation of su(2) and not the full algebra 5u{2F + 1). To create a controllable 
system we need to consider the laser light shift interaction 

^LASER = -^E,{tyEj{t)a,j. (3.7) 



Here, Olj^j IS the polarizability tensor 
A 



a = -J]^=^ (3i 



9,e 



We can reduce this to a more manageable form by expressing the hght shift 
Hamiltonian in terms of its irreducible spherical components 

^...a ^ -\ (.'°'|EP . . ,E- X E, . - 1|EP.,)) . (3.9) 

We can rewrite this Hamiltonian as an effective operator on the atomic spin by 
expressing aij in terms of the generators of angular momentum on F like 

//laser = c(°)|Ep + c(i)(^^) • F + c(2)(|E ■ Fp - \\n'\P?)- (3-10) 

The constants, c^^\ can be found through the Wigner-Eckart as in [59]. 

For our control system we use monochromatic light with polarization along the x- 
direction. In this case the light shift Hamiltonian, dropping constant terms, reduces 
simply to 

//laser = c^^^n^Fl = p^sFl (3.11) 
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Here we can rewrite the constants in terms if the photon scattering rate, 7^, and a 
dimensionless parameter (3 which is a measure of the timescales for coherent versus 
incoherent evolution. Its value depends on the atomic structure and the frequency 
of the driving field and for Cs takes on a maximum value f3 = 8.2 when tuned 
between the hyperfine transitions of the Di line at 894nm. This is enough to allow 
considerable coherent manipulation. Due to technical concerns the laser was an 
"always on" interaction leading to a final control Hamiltonian 

H = + Vt,{t)F, + ^y{t)Fy. (3.12) 

is itself a rank-2 operator of angular momentum, and so we see that this system is 
controllable by direct application of Thm. 1. In the experiment, the photon scattering 
rate is typically around 75/27r 0.77 kHz and the amplitudes of the applied magnetic 
fields are about B ^ 40mG which leads to f^js/"^^ ^0.5 kHz and Q/27r ^ 15 kHz. 

In addition to the nonlinear light-shift, the laser interaction also leads to spon- 
taneous photon scattering. This is important since, in the large detuning limit, the 
photon scattering rate has the same scaling with respect to the intensity and detun- 
ing of the laser as the nonlinear contribution to the Hamiltonian. By choosing the 
optimal parameters we can get some nontrivial evolution before we lose too much 
coherence to spontaneous emission, but since the incoherent and coherent rates are 
intrinsically related there is an upper bound on the length of the coherent control 
fields it is possible to consider. 

3.2 Microwave and rf magnetic fields 

An alternative route to controlling the atomic spins is to employ solely magnetic 
interactions, and remove the necessity of the laser-induced AC-Stark shift. This 
approach has the advantage that we can perform control on the entire electronic 
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ground state rather than one irreducible manifold, a 16-dimensional Hilbert space. 
Additionally, none of the control fields are intrinsically tied to decoherence, with 
spontaneous scattering of rf or microwave photons completely negligible, in principle 
allowing for much richer landscape of possible controls. 

The Hamiltonian describing the atom and its interaction with external magnetic 
fields takes the form given in Eq. 3.1, with laser coupling set to zero. In this control 
scheme we consider the application of three fields, B(t) = Boe^ + Brf(t) + B^w(^)- 
The static bias field Bq defines the quantization axis and Zeeman splittings between 
the magnetic sublevels. The terms Brf(t) and B^w(^) describe magnetic fields os- 
cillating at radio and microwave frequencies, respectively. The hyperfine coupling 
between spins provides an eflFective nonlinearity that will allow full controllability of 
the Hilbert space for appropriate choices of external fields. 

In the linear Zeeman regime, /llbBq ^ A^ the static field acts separately in the two 
irreducible subspaces, and according to the Lande projection theorem, the Hamilto- 
nian is approximately, 

^Mb5Z^/Bo-F(^). (3.13) 

f=± 

Here F^^^ = P±FP± refers to the total angular momentum operator projected onto 
the subspaces with quantum number F±. Neglecting the nuclear magneton contribu- 
tion, the g-factors for the two manifolds have equal magnitude but opposite sign, i.e. 

= —g^ = 1/F+. The hyperfine coupling plus bias magnetic field thus determine 
the static Hamiltonian, 

Ho = (P+ - P-) + r2o(Fi+) - Fi")), (3.14) 

where AE^f = AF^ is the hyperfine splitting and = /llbBq/F^ is the Zeeman 
splitting between neighboring magnetic sublevels. 

As our first control field, we consider rf-magnetic fields oscillating near the fre- 
quency of the Zeeman splitting, cjrf ^ ^^o, realized by Helmholtz coils driven with 
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Figure 3.2: The ground state hyperfine manifold of ^^^Cs. Rf-magnetic fields (in 
red) lead to independent rotations in the two manifolds. Microwaves (in blue) are 
the generators of rotation in a two-dimensional subspace between states in the two 
manifolds, here the stretched state transition |4,4) |3,3). 

the appropriate current. We take two sets of coils that produce fields with x and y 
polarization, independent amplitude and phase control, but equal carrier frequency, 
cjrf- Again, for a moderate current such that the amplitude of the magnetic field 
is in the linear Zeeman regime, the rf-Hamiltonian takes a form equivalent to the 
interaction with the static field 



The time dependent amplitudes {Qx{t)^Qy{t)) and phases (^^^(t), (/)^(t)) of the two 
sets of rf coils will be used to control the system. 

To better understand the eflFect of the rf field, consider a resonant interaction, 
cjrf = ^0- In the rotating frame, H^i{t) H^^{t) = Ul^H^i{t)U^i^ where U^i = 



= ^x{t) COS {uo^it 

+ ^y{t) COS {Uj^ft 



(3.15) 
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z-axis m 



exp |— i6Jrft(i^i^^ — is a rotation of the two manifolds about the 

opposite directions, fJ^^ fJ^^ cos(6Jrft) ± F^^^ sin(cL;rft), F^^^ Fy^^ cos{ujj.{t) =F 
fJ^^ sin(cL;rft). Performing this unitary transformation and averaging over a cycle, 
the rf-Hamiltonian in the rotating wave approximation is. 



H^At) - ^cos(0.(t))(FW-Fi-)) 
^sin(0.(t))(FW+F(-)) 
^cos(0,(t))(FW-F(-)) 
- ^sin(0,(t))(FW+Fi-)). (3.16) 



2 

n 

+ 2 

n 

+ 



Rf-control of the two spin manifolds differs from the familiar spin resonance problem. 
In the latter, a single magnetic field in either the x or ^/-direction would be sufficient 
to generate the entire SU (2) algebra for rotations. With two irreducible manifolds 
there is an added freedom - the two angular momenta F+ and F_ can rotate in the 
same or opposite directions. Amplitude and phase control of two rf-magnetic field 
polarizations allows us to perform arbitrary and independent rotations on the two 
hyperfine manifolds. With only a single direction of Brf we would be restricted to 
either co-rotating or counter-rotating in the two subspaces. 

The weak rf-magnetic fields alone will not be sufficient to fully control our atomic 
system; they don't couple the F+ and F_ manifolds, nor do they provide a nonlinear 
Hamiltonian within these subspaces. In order to make our system fully controllable, 
we look to resonant microwaves. While the fundamental Hamiltonian governing the 
microwaves is exactly of the same form as the quasistatic magnetic fields, the resonant 
behavior leads to very different dynamics than the previous interactions. Depending 
on the polarization and frequency, the microwave couples a Zeeman sublevel in F+ 
manifold with one in the F_ manifold whose magnetic quantum number differs by 
Am = 0, ±1. For a sufficiently strong bias Bq we can ignore any off-resonant exci- 
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tation, and restrict the Hamiltonian to act only on a 2D subspace spanned by the 
states we are trying to couple. In that case the microwave Hamiltonian has the form 



Hj^^{t) = ri^wW cos {u^^t - (t)^^{t))a^, (3.17) 

where is the Pauli sigma-x matrix for this pseudospin, cr^ = m+)(F_, m_| + 
m_)(F+, m+l and Vt^^{t) is the (time-dependent) Rabi frequency depending on 
the microwave power and the transition matrix element. Again, the amplitude and 
phase of the microwave fields are control parameters. In this subspace, the problem 
takes the form of the standard two-level resonance problem. We must take care 
in going to the rotating frame to account for the simultaneous transformation we 
perform due to the rf-fields. The complete frame transformation is achieved by the 
unitary matrix 

U = f/rf exp (P+ - P-) 1 , (3.18) 



2 

where a — uo^^ — (m+ + m_)cc;rf. Under this transformation, the Hamiltonian in the 
rotating wave approximation for resonant microwaves is 

+ ^^sm{4>,Ut)W (3.19) 

generating rotations of this pseudo-spin on the Bloch sphere. 

Combining the static, rf, and microwave interactions the final Hamiltonian in the 
rotating frame is 

H'{t)^H'^ + H',,{t) + H'^^{t). (3.20) 

Allowing for a finite detuning of the oscillating fields from resonance, the static 
Hamiltonian in the rotating frame becomes, 

^0 = % {P+ - P-) + Arf(Fi+) - Fi")), (3.21) 
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where A^w = ^/iw — ^Ehf — {m + m')ijj^i is the effective detuning of the microwaves 
from the two-level transition of interest, |F_,m_) |F+,m+), and Arf = uOri — 
is the rf detuning. This, together with Eqs. (3.16,3.19), defines the Hamiltonian we 
employ for control, and which we will analyze for use in arbitrary state preparation. 



For this Hamiltonian system, with arbitrary control of the amplitude and phase 
of the two orthogonal sets of rf-coils and a single microwave field, the control algebra 
generated by the six operators {Fx'^\ Fy'^\ Fx~\ Fy~\ (Jx^ (Jy) is 5u(rf) in its entirety. 
In this case, it is possible to prove controllability analytically for an arbitrary alkali, 
with an arbitrary nuclear spin /. 



The proof is as follows, first we would like to show that with our Hamiltonian the 
subspaces and F_ are independently controllable. To show controllability of the 
manifold we require an operator that has a nonzero overlap with a rank-2 tensor 
on that space. Restricted to the subspace, the operator looks like a projector 
onto some particular sublevel, |F+, m+)(F+, m+|. The overlap of this projector with 
T^^ is Tr ^|F+, m+)(F+, m+|To^^^ = ^5/11 m+|2, 0; F+, m+), which is nonzero 
for all values of m+. Of course, also has support in the F_ manifold, however, 
Fx'^\(jz does not. Since commuting by Fi^^ can't change the rank of a tensor, 
we are left with an operator confined to the manifold that has a nonzero overlap 
with some rank-2 tensor, and so according to theorem 1, we have complete control 
of the manifold. This proof directly carries over to the F_ manifold. 



At this point we have shown that we have full controllability over both the 
and the F_ subspaces, as well as the 2-dimensional subspace coupled by the resonant 
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microwaves. We can write this in matrix form 

/I I \ 



v 









S2 



(3.22) 



where we have ordered the basis vectors so that the states coupled by the mi- 
crowaves are adjacent to each other. We have shown that we can simulate any 
operator that only has matrix elements within the three boxes in Eq. 3.22, i.e. all 
operators that have support only on the diagonal and super diagonal matrix ele- 
ments. The irreducible representations of angular momentum, Jx and J^, on the 
entire space have support only on the super-diagonals. Therefore, we can simulate 
Jx and Jy. According to theorem 1 all we need to show for controllability is that 
we can simulate some operator with a nonzero overlap with some rank 2 operator. 

(2) 

Since we can simulate any diagonal operator, we can simulate Tq . It thus follows 
that {F^'^\Fi^\ F^~\f^~\(Jx, (Jy) generates su(rf). 

Though sufficient, the entire available set is not necessary to achieve controllabil- 
ity. In practice, one can reduce the number generators in the control algebra and still 
implement an arbitrary unitary. For an experiment, it is important to understand 
which components are really necessary so that we can chart the tradeoffs between 
ease of implementation and controllability. In order to study the capability of various 
reduced sets of controls, we resort to numerical approach discussed in Ch. 2.1. 

We carried out this procedure for the specific example of ^^^Cs with nuclear spin 
/ = 7/2 to study the capability of a variety of control sets to generate the entire 5u(16) 
algebra. We considered 8 different microwave configurations: controlling or fixing 
the amplitude and the phase of the fields, and whether or not we are detuned from 
resonance. The two cases where both the amplitude and phase are controlled and 
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1 fiw trans, 
fixed ampl, 
fixed phase, 
resonant 


1 /xw trans, 
controlled ampl, 
fixed phase, 
resonant 


1 fiw trans, 
fixed ampl, 
controlled phase, 
resonant 


1 /xw trans, 
fixed ampl, 
fixed phase, 
detuned 


1 fiw trans, 
controlled ampl, 
fixed phase, 
detuned 


1 /xw trans, 
fixed ampl, 
controlled phase, 
detuned 


1 rf polarization, fixed ampl, 
fixed phase, resonant 


# 


# 


O 


# 


O 


o 


1 rf polarization, controlled ampl, 
fixed phase, resonant 


W 


W 


o 


o 


o 


o 


1 rf polarization, fixed ampl, 
controlled phase, resonant 


1 — 1 


1 — 1 


w 


w 


w 


w 


1 rf polarization, fixed ampl, 
fixed phase, detuned 


# 


1 — 1 


w 


W 


w 


w 


1 rf polarization, controlled ampl, 
fixed phase, detuned 


1 — 1 


1 — 1 


9 


w 


CJ 


U 


1 rf polarization, fixed ampl, 
controlled phase, detuned 


1 — 1 


1 — 1 


w 


w 


9 


w 


2 rf polarizations, fixed ampl, 
fixed phase, resonant 


# 


# 


o 


# 


o 


o 


2 rf polarizations, controlled ampl, 
fixed phase, resonant 


□ 


□ 


• 


• 


• 


• 


2 rf polarizations, fixed ampl, 
controlled phase, resonant 


□ 


□ 


• 


• 


• 


• 


2 rf polarizations, fixed ampl, 
fixed phase, detuned 


• 


□ 


• 


• 


• 


• 


2 rf polarizations, controlled ampl, 
fixed phase, detuned 


□ 


□ 


• 


• 


o 


o 


2 rf polarizations, fixed ampl, 
controlled phase, detuned 


□ 


□ 


• 


• 


o 


o 



Table 3.1: Table exploring controllability of the system for a variety of configura- 
tions: one microwave field driven on different two-level transitions, |F = 3, M) 
\F = A^M')^ amplitude and/or phase control, one or two sets of orthogonal rf 
coils (rf polarizations), and resonant vs. detuned fields. The different configura- 
tions yield one of four different outcomes: (green circle) all microwave transitions 
provide full controllability, (yellow square) all transitions but the clock transition 
|3,0) |4,0) provide full controllability, (orange pentagon) only the transitions of 
the form |3, ±3) |4, ±4) and |3, ±3) |4, ±2) provide full controllability, and (red 
octagon) no transitions yield controllable Hamiltonian dynamics. In this calculation 
we consider all valid microwave transitions, not only the energy non-degenerate ones. 



where the amplitude is fixed but the phase is controlled can be shown to be equivalent. 
In the rf configurations we also allow for one or two orthogonal sets of magnetic coils. 
The last free parameter is the choice of which microwave transition we excite. We 
assume arbitrary frequency and polarization selectivity of the desired transition for 
this purpose. The results are summarized in summarized in table 3.1. In each box we 
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enumerate the set of microwave transitions that yield controllable dynamics. We find 
that our system is controllable in a wide number of regimes, though there are some 
configurations in which it is not. For example, out of all the choices for microwave 
transitions, the clock transition, 0) 0), is controllable in the least number 

of scenarios. This shouldn't come as much of a surprise since we are controlling the 
system with rf magnetic fields and this transition's insensitivity to magnetic fields is 
what makes it useful for precision metrology. 

It is interesting to note that there exist configurations that are controllable in 
which there is one time-dependent control waveform and some fixed time-independent 
interactions. This is the simplest system one could expect to find, and allows for 
bang-bang control, a well-studied protocol. In the next chatper, however, we look 
at the control systems that utilize more parameters, decreasing the time needed for 
state preparation. 
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Chapter 4 

State Preparation with AlkaH 
Atoms 



In this chapter I will discuss the application of the theoretical methods for state 
preparation discussed in Ch. 2 to the the physical systems of alkali atomic spins 
discussed in Ch. 3. In the context of control via AC-Stark and quasi-static magnetic 
fields, this protocol was carried out in the laboratory and shown to yield good results 
as described below. In the context of microwave/rf control, I have devised new 
protocols for control which have been studied numerically. Experimental test should 
be forthcoming in the near future. 



4.1 A state preparation algorithm 

We seek to design Hamiltonian evolutions that take an initial known pure quantum 
state to an arbitrary pure state in the Hilbert space of interest. We would like to 
maximize fidelity as a functional of the control waveform given by 

F[c{t)] = \{iJtarget\U[cit)]\^o)\'. (4.1) 
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As we explored in chapter 2.2, this problem has an extremely favorable topology 
devoid of local maxima, and therefore, a local search of the space of control fields, 
starting from any random initial guess, will find a global maximum of the fidelity. 
For this problem, gradient searches perform about as well as more computationally 
intensive searches like genetic or simulated annealing algorithms. 

In a real system, we will violate some of the assumptions required for the proof 
that there are no local maxima. There will always be some decoherence and one 
does not have infinite time to perform the control. In fact, we would like to perform 
state preparation as fast as possible in order to combat decoherence and various 
inhomogeneities that lead to accumulated errors. Additionally, we need to consider 
control fields that have a limited bandwidth and slew rate constraints. For these 
realistic conditions, not every gradient search from an arbitrary starting point yields 
a global maxima. Nonetheless, we have found empirically that the results of the 
theorem are approximately true with moderate decoherence and after a sufficient 
time. We still find excellent protocols after making only a small handful of searches, 
and these can be further filtered to find control waveforms that perform well under 
realistic operating conditions. 

As we are dealing with the optimization of waveforms that are functions of contin- 
uous time, the first step is to transform the problem into a search for a finite number 
of values at discrete times. The physical constraints of bandwidths and slew rates of 
the controllers provide a natural scale. There is a minimum interval during which a 
field can vary over a maximum range. A discretized version of a control waveform 
is thus specified as a vector of values within this range at these fixed intervals. The 
continuous control waveforms are then found by interpolation using cubic splines, 
consistent with the bandwidth constraints, at least on a fine enough grid for use in 
our numerical integration of the Schrodinger equation. 

We create optimal control waveforms by first fixing the total time of the state 
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preparation procedure. Due to our discretization technique, fixing the total time 
fixes the number of optimization variables. Starting from a randomly chosen initial 
vector of control waveform values, bo, we perform a gradient ascent search by taking 
small steps in the direction of steepest ascent, i.e. 



An optimal value corresponds to the maximum, where the gradient approaches zero. 
We performed this search numerically on a Matlab cluster by optimizing waveforms 
from a handful of random seeds in parallel, and then chose the one that gave the 
highest fidelity. The actual gradient search itself was performed using a canned algo- 
rithm from Matlab's "Optimization Toolbox." An alternative approach would have 
been to use the Gradient Ascent Pulse Engineering (GRAPE) algorithm developed 
in [17]. While this algorithm has been used to great affect in a number of quan- 
tum control protocols, it was of no use in our scheme for controlling atomic spins. 
For completeness, I give a brief summary of the GRAPE algorithm here and its 
limitations. 

The GRAPE algorithm is a gradient search algorithm whose outstanding feature 
is that the gradient is computed in a way that is much more efficient than a stan- 
dard numerical differentiation routine. In the simplest incarnation of the GRAPE 
algorithm we imagine the controls c to describe the amplitude of square pulses with 
one time-varying control field. We can write the fidelity of state preparation as 



When we numerically evaluate J we need to perform matrix multiplications. If we 
were to try to we calculate the A^-dimensional gradient numerically we would need to 
evaluate J, approximately times. This leads to a scaling of O(A^) for computing 
the gradient numerically. 

The GRAPE algorithm allows us to compute the gradient with a cost of 0{N). 



(4.2) 




-iAt(Ho+CNHi) -iAt(Ho+CN-iHi) -iAt(Ho+ciHi) \-\\2 



(4.3) 



54 



Chapter 4. State Preparation with Alkah Atoms 



The trick is to simultaneously forward evolve the initial state while backwards evolv- 
ing the final state. We define two sets of vectors 

= ^-iAt{Ho+CjHi) ^-iAt{Ho+Cj-iHi) ^ ^ ^ ^-iAt{Ho+ciHi) ^ j^^^^^j 

with |i(0)) = \i) and (/(r)| = (/|. Clearly 

J = Umimi' = K/(j)|e-^^'^''°^^^''^^N(j - I))!', V < J < T. (4.6) 

In the GRAPE algorithm we first compute (/(1)| and at a cost of 0{N) matrix 

multiplications. If At is sufficiently small we can compute 

J -iAt(Ho+ciHi) 

= (4.7) 

at a cost that is constant in A^. To compute dJ/dc2 we require the vectors (/(2)| 
and \i{2)) which we obtain by evolving 

(/(2)| = (/(l)|e'^*(^°+^^^^\ |i(2)) = e-'^*(^°+^^^^)|z(l)). (4.8) 

This only takes two matrix multiplications. We can repeat this N times in order to 
get all the components of VcJ with only 0{N) matrix multiplications. 

The assumption I have made in this simple description of the GRAPE algorithm 
is that the number of optimization variables corresponds to the ratio between the 
total time of the state preparation and the sampling time for the Schrodinger inte- 
grator. This assumption is valid in the case where the maximal slew rates for the 
applied fields are much larger than the Rabi or Larmor frequencies of the fields, as 
is the case in liquid state NMR systems. In more concrete terms, the restriction 
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for convergence in the GRAPE algorithm is that the integration step length of the 
Schrodinger integrator At must satisfy 

At<^\\H{t)\\-' (4.9) 

where, is the seminorm of H{t)^ or the largest frequency component of 

the Hamiltonian. This constraint arises from the requirement that the derivative 
of exp (— zAt(i/o + CiHi) is approximately (— zi/i)exp (— zAt(i/o + CiHi)). For our 
state preparation routine the constraint is that 

At'«||^||-\ (4.10) 

requiring that the Hamiltonian must be relatively uniform over a time step. This 
leads to a computational cost of 0{T/At) for the GRAPE algorithm and a cost of 
0{{T / At'Y) for ours. For instances of the control systems in this dissertation it 
turns out that the second term is smaller and thus the GRAPE algorithm is not an 
efficient approach for our system.. 

For the experimental implementation, we additionally require that the state 
preparation protocol is at least somewhat robust to inhomogeneities and noise. To 
enforce this we use a two round optimization. First, we find a set of state preparation 
protocols using the above technique. For the parameters we considered this would 
typically yield fidelities of greater than 0.99. At this point we switch to a more re- 
alistic estimate of control performance by modeling the evolution with a full master 
equation that incorporates decoherence from light scattering and inhomogeneity of 
the nonlinear strength across the atomic ensemble. This allows a second stage of 
optimization starting from the waveform generated in round one and using the more 
complete but computationally intensive model to predict the yield, which is now de- 



fined in terms of the overlap y — Try P^P' PtPp^'^ between the target density matrix 
prj. and the predicted density matrix pp. 
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4.2 Quasi-static magnetic fields and light shift 



We demonstrate, in [ ], the quantum control of the spin-angular momentum associ- 
ated with the F = 3 hyperfine ground state of individual ^^^Cs atoms, a 2F + 1 = 7 
dimensional Hilbert space. Starting from an easily prepared fiducial state we use 
time-dependent magnetic fields and a fixed AC Stark shift to design and implement 
near-optimal controls and produce a range of target states. We evaluate our con- 
trol performance by experimentally reconstructing the entire spin density matrix [43] 
and computing the overlap between the measured and target states. In most cases 
the estimated yield is in the 0.8 — 0.9 range, limited by errors in the control fields 
and to a lesser extent by decoherence from light scattering. The measured states 
can be compared also to the predictions of a full model that includes the effects of 
errors and decoherence. Typical fidelities between measured and predicted states 
are around 0.9, which is close to the resolution limit of our procedure for quantum 
state estimation. We further use this universal approach to generate spin-squeezed 
states and compare against a method based on adiabatic evolution [ ]. The latter 
is more robust against errors in the control fields, but also slower and thus more 
sensitive to light scattering and decoherence. Large spins provide a testing ground 
for the design of accurate and robust controls in a system where the Hamiltonian is 
well known and where errors and dissipation are well understood and can be accu- 
rately modeled. From a practical perspective, quantum control of hyperfine states 
has direct relevance to proposals for neutral atom quantum computing [ ] wherein 
qubits (or higher dimensional qudits [ ]) are encoded in the ground-state manifold, 
and may provide a simple route to modest spin squeezing and accompanying gains 
in precision atomic magnetometry [62]. 

The combination of a time-dependent magnetic field and a constant x-polarized 
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Figure 4.1: Quantum control of a large atomic spin, (a) Schematic of the experiment, 
(b) Example of a control waveform (/>(t). (1-4) Wigner functions at four stages during 
the control sequence. Both sides of the sphere are shown. The final result is close to 
the target state \iptarget) = = 2) + |m^ = — 2))/\/2. (c) Density matrix (absolute 
values) and Winger function for \iptarget) 

light field, discussed in chapter 3.1, results in a control Hamiltonian [ ], 

Hc{t) = gffiBB{tyF + /?7.Fi (4.11) 

A schematic of our setup for spin quantum control is shown in Fig. 4.1(a). We 
begin with a sample of a few million Cs atoms, captured and laser cooled to ^ 
in a magneto-optical trap and optical molasses. Once the atoms are released from 
the optical molasses their spin state is initialized by optical pumping into a state 
of maximum projection along the ^-axis, {ipo) = \F = S^rrty = 3). We drive the 
spins by applying a time-dependent magnetic field from a set of low-inductance coils 
driven by arbitrary waveform generators, and by applying a static light shift from 
an optical probe beam. Using an all-glass vacuum cell, avoiding nearby conductive 
or magnetizable materials, and synchronizing our ^ 0.5ms duration experiment to a 
fixed point during the AC line cycle allows us to null the background magnetic field 
to a few tens of /xGauss without the use of shielding or active compensation. The 
applied magnetic field can be controlled with an accuracy better than one percent in a 
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bandwidth of more than 100 kHz. Immediately foUowing a period of quantum control 
we estimate the resulting quantum state as described in [ ] . In this procedure the 
control magnetic and optical fields are applied to drive the spins for an additional 1.5 
ms, while continually and weakly measuring a spin observable through its coupling 
to the probe polarization. To reduce the effect of noise, the measurement signal is 
averaged over 16 repetitions of the experiment and the full density matrix determined 
from the measurement record and the known evolution. 

The objective is to start from the state lipo) and to produce a specified target state 
li^target) by modulatiug the field B(t) for a fixed time r. With readily available mag- 
netic fields the timescale for geometric rotations is much shorter than for nonlinear 
evolution driven by the light shift, and the latter therefore becomes the time-limiting 
element of most transformations. In our experiment the maximum available Larmor 
frequency is 15 kHz and the nonlinear strength is P'js ^ 27r x 500 Hz. Under these 
conditions there is no significant sacrifice in control performance when the set of 
available rotations is somewhat restricted. We therefore choose the magnetic field 
to have constant magnitude and time-dependent direction in the x-y plane. With 
this simplification the control Hamiltonian is completely determined by the time 
dependent angle (/){t) between B(t) and the x-axis. 

B{t) = Bi (cos e, + sin e^) (4.12) 

The state \iptarget) and the transformation {ipo) I'lptarget) belong to a = 7 di- 
mensional Hilbert space and can be specified by a set of 2rf — 2 = 12 real numbers, 
and full control therefore requires at least that many free parameters in the control 
Hamiltonian. To ensure sufficient ffexibility we specify the control waveform (j){t) by 
its values {0^} at = 30 discrete time steps. 

An example of an optimized control waveform is shown in Fig. 4.1(b), along with 
Wigner function representations of the spin wavepacket [ ] (see appendix A) at a few 
steps during the transformation as calculated using the complete master equation. 
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Figure 4.2: Examples of target and measured density matrices (absolute values). 
The target states are (a) (Im^ = 2) + \mz = — 2))/\/2 , {h)\mx = 2), (c) T,ymy\my) . 
The experimental yield is indicated for each case. 

Note that the nonlinear evolution initially produces a squeezing ellipsoid which later 
wraps around the sphere so that interference effects can be manipulated to create the 
desired state. The end product is very close to the target state shown in Fig. 4.1(c). 
According to our model this and a wide variety of other control waveforms all produce 
yields near 0.95. Taking into account imperfect optical pumping in our experiment 
(the initial population in l^po) is ~ 0.96) reduces the expected yields to around 0.90. 

We have generated and tested a sample of control waveforms designed to produce 
21 different pure spin states. Fig. 4.2 shows three examples of target and measured 
density matrices, with yields falling in the range 0.87-0.97. A more complete statis- 
tics of yields for over a hundred experimental realizations of control is compiled in 
the form of a histogram in Fig. 4.3(a), showing a fairly broad distribution centered 
on respectable value of 0.8. It is also informative to compare the experimentally mea- 
sured density matrices against the density matrices pp predicted by our model, 
as quantified by the fidelity T = Tr^ Pp^^PuPp^^- ^ig- 4.3(b) shows a histogram 
of fidelities for our data set. Note that both yield and fidelity can be affected by 
control errors (the real state is different from pp) as well as state estimation errors 
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Figure 4.3: Histograms of (a) yields and (b) fidelities of measured vs. predicted 
states, (c) & (d) Yields and fidelities when each measured state is geometrically 
rotated to optimize overlap with the predicted state. 



(the real state is different from p^), and that there is no way to distinguish between 
these possibilities. Numerical modeling shows that small background magnetic fields 
or miscalibration of the control fields will lead to apparent geometric rotations of 
the final state, but such errors are too small in our experiment to significantly affect 
the outcome. The obvious outliers in the yield and fidelity distributions are asso- 
ciated with two specific control waveforms, and closer examination shows that the 
estimated states are rotated relative to the predicted states. The axis of rotation 
corresponds the direction of the magnetic field at the transition between the control 
and state estimation phases, which suggests a problem with the way the correspond- 
ing control waveforms were joined together. We can numerically rotate a given 
to maximize its fidelity relative to pp and obtain new values for yield and fidelity. 
Carrying out this procedure for all data points takes care of the outliers without 
otherwise changing the yield distribution significantly, as shown in Fig. 4.3(c). This 
distribution can reasonably be interpreted as a measure of our ability to control the 
spins in a well designed experiment. The fidelity distribution (Fig. 4.3(d)) remains 
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Figure 4.4: Spin squeezing by adiabatic control, (a) Normalized squeezing parameter 
vs. final magnetic field for the squeezed and anti-squeezed components. Dashed 
lines: perfect squeezing. Open symbols: predictions from our theoretical model. 
Filled symbols: experimental results, (b) Target and (c) measured Wigner functions 
corresponding to the smallest observed ^. 

peaked at 0.9, which we know from experience to refiect the accuracy of our state 
estimation algorithm. Finally we note that random errors in state estimation are 
far more likely to decrease than increase the apparent yield. A simple error model 
based on Gaussian random displacements in state space indicate that the yields are 
probably 10% larger on average than indicated by Fig. 4.3(d). This puts most yields 
in the range 0.8-0.9, in good agreement with the ~ 0.9 predicted by the model used 
to design the control waveforms in the first place. 

To further explore quantum control in our system we have studied the generation 
of spin squeezing both by optimal control as outlined above and by the adiabatic 
scheme described in [ ]. The latter begins with an initial state, {ipo) = |F = 
S^niy = —3), which has equal uncertainties for the components AFx and AF^ and 
is often referred to as a spin-coherent state. This state is a good approximation to 
the ground state of the control Hamiltonian Hc{t) when the magnetic field is of the 
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form B(t) = B{t)y and B{t) is large. As the field magnitude is slowly reduced the 
state adiabatically evolves so as to minimize the squeezing parameter ^ = AFx/\{Fy) \ 
of relevance for metrology [ ]. Fig. 4.4(a) shows the progression of squeezing and 
anti-squeezing relative to a spin-coherent state with the same \{Fy)\. Up to ^ 4 dB 
of squeezing is seen in the experiment, in good agreement with the predictions of our 
model. For the small spin magnitude used here the squeezing is quickly limited by the 
decrease in \{Fy)\ as the squeezing ellipse wraps around the sphere. Fig. 4.4(b)-(c) 
shows Wigner functions for the target and actual state for the smallest ^ achieved 
in our experiment (~ 80% of the coherent state value). We have produced the 
same spin squeezed states via optimal control, with small but significant reductions 
in both squeezing and yield. This suggests that gains from reduced decoherence 
(optimal control is as much as five times faster) is oflFset by increased sensitivity to 
control errors. 



4.3 Microwave and rf magnetic fields 



In [18], we developed the microwave and rf control system in chapter 3.2, and ap- 
plied our state preparation technique to the complete 16-dimensional ground state 
manifold of ^^^Cs. We take a static bias field to produce a Zeeman splitting of 
= 1.0 MHz, sufficient to give excellent resolution of the magnetic sublevels, 
but well within the linear Zeeman regime. The rf field power is chosen so that 
on resonance the rotation rate is characterized by Qj.^ = 15 kHz. As a generic 
case, we take one microwave field, resonant on one of the stretched transitions 
|F = 3,M = ± — 3) |F = 4, M = ± — 4), where the microwave Rabi fre- 
quency is largest, and the system is controllable in a wide variety of scenarios. The 
microwave power is chosen to give a Rabi frequency f^^w = 40 kHz. The slew rates 
constrain the maximum rate of change of amplitude and phase of the control fields. 
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In the case of the rf-magnetic field, a "slew time" of Trf = 10/xs fixes the slew rates 
on the amplitude to 1.5 kHz//is and phase to 0.2 7r//xs. In the case of microwaves, 
faster control is possible, with a slew time of r^w = 1.0 //s, or amplitude and phase 
slew rates of 40 kHz//xs and 2.0 7r//xs respectively. 

Two examples of the end product of this optimization are shown in Figs. (4.5,4.6) 
for target states :^(|4, 4) + 13, -3)) and ;^|4, 4) + 1(|3, 3) + |3, -3)) respectively. The 
initial state for these examples is the stretched state |4, 4), a state easily reached by 
optical pumping. We control the amplitudes and phases of rf coils in both the x and 
y directions, as well as the amplitude and phase of a resonant microwave that couples 
the states |4, —4) and |3, —3). In Figs. (4.5,4.6) we show the Cartesian components 
of the three control fields cos (j) and Q sin 0) over the entire state preparation time 
of 150/xs. The figures show snapshots of the evolved state at five different times, 
identified as times (0)-(4). Two different representations of the state are shown: bar 
charts of the absolute values of the density matrix elements, and a generalized spher- 
ical Wigner function. The spheres on the diagonal represent the Wigner functions in 
the irreducible subspaces and the off-diagonal spheres represent the coherences 
between the manifolds. For details see Appendix A. The fidelities of preparation in 
both cases are greater than 99%. With a state preparations time of 150/xs moderate 
searches yield high-fidelity waveforms. More intensive optimizations can yield faster 
control waveforms. 

Our gradient search algorithm leads to waveforms that cause the system to un- 
dergo quite complex dynamics, as evidenced by the intermediate states seen in the 
course of the evolutions. Figs. (4.5,4.6). One may wonder whether there are simpler 
choices, since given a fixed initial state, there are many different waveforms that 
lead to same target state. While our method does lead to waveforms that are hard 
to intuitively understand, some recent studies [ ] suggest that the waveforms de- 
rived from gradient searches may be more robust than those that come from more 
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Figure 4.5: ^{\4^A) + |3, -3)) prepared with fidelity 0.993. 
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Figure 4.6: ^|4,4) + |(|3,3) + |3, -3)) prepared with fidelity 0.995. 
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geometric algorithms. 

We discussed the mathematical conditions necessary for our Hamiltonian dynam- 
ics to be controllable. These conditions, while useful for ruling out large classes of 
Hamiltonians as unsuitable for our purposes, tell us nothing about the relative per- 
formance of different control scenarios. Our figure of merit is the time after which 
we can be reasonably sure that our optimization will find a high fidelity waveform 
for any target state. To determine this time for a given control protocol, we run 
our optimization up to a given final control time over a large collection of randomly 
chosen states and determine the average fidelity. In this section we examine these 
results and discuss some of the tradeoflFs and bottlenecks that might be encountered 
in the lab. 

There are many parameters in this system that we can manipulate, including 
the number of independently controlled rf polarizations, the number of resonant mi- 
crowave frequencies, the types of controls (amplitude vs. phase), detuning, slew 
rates, and the strengths of the diflFerent fields. Based on some of our previous experi- 
ments we set as a baseline: one microwave frequency, two orthogonal rf polarizations, 
rf power giving f^rf =15 kHz, a microwave Rabi frequency of f^^w = 40 kHz, a rf slew 
time of 10 //s, and a microwave slew time of 1.0 /xs. While we could independently 
vary all these parameters, this would be an unwieldy computation. Here we fix some 
of the parameters that are unlikely to diflFer in the future experiments we are con- 
sidering. In particular, we fix the rf slew time to be 10 /xs and consider control with 
two sets of rf coils. For simplicity we also consider all fields to be resonant, and the 
microwaves to couple the stretched states. 

Statistics were collected by running the state preparation algorithm for 10 differ- 
ent random states found by sampling using the Harr measure on SU (16) [ ]. In all 
cases the initial state was the |4, 4) state. For each combination of total time, target 
state, and system configuration, we run the optimization 20 times starting from dif- 
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Amp vs. Phase Control; MW Slew Times Larmor and Rabi 




Figure 4.7: Plots of the average fidelity of state preparation for different control 
configurations and total preparation times. Each point represents the fidelity aver- 
aged over a set of 10 states randomly chosen from the Harr measure. For each state 
and configuration, the gradient search was performed with 20 random seeds and we 
chose the protocol that generated the highest fidelity. 



ferent random seeds of the vector that defines the control waveform. Out of this set 
of 20, we choose the highest fidelity preparation. The fidelities from the 10 random 
states are averaged to produce the data points shown in Fig. 4.7. In principle, more 
iterations would yield higher fidelity waveforms, but it is useful to understand which 
types of high-fidelity controls can be found after only modest searches. 

In Fig. 4.7a, we study the effect of varying the characteristics of the microwave 
field. We compare the performance of one vs. two resonant microwave frequencies 
on one or both of the stretched transitions, |3, 3) |4, 4) and |3, —3) |4, —4). In 
addition, we examine the effect of removing control of the microwave amplitude (a 
scenario that still allows for full controllability of the system, as discussed in Ch. 3.2). 
As expected, since the microwave Rabi frequency is larger than the rf Larmor fre- 
quency, increasing the number of microwave fields has a large effect. On the other 
hand, it was surprising that fixing the microwave amplitude, thereby substantially 
decreasing the number of control parameters, yielded higher fidelity waveforms. We 



68 



Chapter 4. State Preparation with Alkah Atoms 

suspect that while there most hkely exist higher fidehty waveforms with control of 
both amplitude and phase, increasing the number of microwave control parameters 
rapidly increases the dimension of the search space, requiring many more iterations of 
our algorithm to find a superior waveforms, on average. This suspicion is reinforced 
by Fig. 4.7b, where we consider the effect of microwave slew time. With our baseline 
parameters, it would appear that increasing the microwave slew time doesn't really 
limit the optimized control performance. In fact, the smallest slew time we consid- 
ered, 1.0 /xs, performed slightly worse than the other slew times, including the case 
where the microwave amplitudes are fixed. As a reminder, the slew time determines 
the information content of our waveforms, and thus the number of optimization vari- 
ables. Again, we see that for the modest searches we are performing, decreasing the 
dimension of the search space counterbalances the loss of control. 

In Fig. 4.7c, we study the effect of the power in the rf and microwave fields. For 
these simulations we fixed the amplitudes of the fields and solely control their phases. 
We find that varying the microwave power around our baseline makes little difference. 
The rf power is slightly more important, but increasing the Larmor frequency above 
the baseline has a fairly small effect. These results indicate that the slew rate and 
bandwidth constraints we have imposed on the rf magnetic fields are the bottleneck 
for controlling the system, and limit the ability to more rapidly control the system 
through increases in power. It would appear that the microwave parameters we 
employ as our baseline are also well above the limits imposed by this bottleneck 
and we can safely reduce the microwave power and slew rates without sacrificing 
performance. The rf Larmor frequency we employ is commensurate with the slew 
rate constraint. 

By optimizing many state preparations for a variety of control configurations we 
find state preparation protocols with this system that take between 50 — 150//S. We 
can compare this to the types of control waveforms that were implemented in our pre- 
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vious work that employed a nonlinear AC-Stark shift to achieve controllability [41]. 
The waveforms we find here are about an order of magnitude faster, control a Hilbert 
space that is double the dimension, and have negligible decoherence as compared to 
the intrinsic decoherence that arises from spontaneous emission. 
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Efficiently Constructing Arbitrary 
Unitary Maps on Qudits 



In chapter 2.2,1 discussed the difference between the problems of unitary construction 
and state preparation. The stochastic search, state preparation techniques used in 
the last chapter are ill-suited to the problem of designing general unitary maps. Most 
known techniques for constructing arbitrary unitary maps fall under the category of 
geometric constructions (Ch. 2.3) which, while powerful, lack broad applicability. 
In [45], my coauthors and I developed a new type of unitary construction protocol 
which is a hybrid of stochastic/geometric construction, similar to the protocol in [^^^]. 
Essentially, we leverage off of our ability to efficiently generate state preparations, 
and then splice state preparations together in a geometric way to create a general 
unitary map. The types of Hamiltonian dynamics that this construction applies to 
have some restrictions beyond controllability. These restrictions are, however, much 
less stringent than those in most geometric techniques. 

In chapter 5.1 I present our hybrid protocol for constructing general unitary maps 
by combining efficient numerical searches with a deterministic algorithm. In addition 
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to unitary maps on the full Hilbert space, this scheme allows us to construct maps 
on a subspace with a complexity that scales as the dimension of that space. In 
chapter 5.2, our unitary matrix construction is applied to control the large manifold 
of magnetic sublevels in the ground electric states of an alkali atom (e.g. -^^^Cs) 
[18]. We show how to construct a set of unitary matrices on SU{d) that are often 
considered as qudit logic gates in a fault-tolerant protocol. In addition, we apply 
our construction for subspace mapping to encode logical qubits in our qudit, and 
simulate an error correcting code that protects against magnetic field fiuctuations. 



5.1 Unitary construction 

In this section we define an efficient protocol for constructing arbitrary unitary maps 
based on state preparation. Any unitary matrix has an eigen-decomposition, 

/7 = 5^e-^^^>,)((/),| =[]e-^^^-l^^-><^^-|, (5.1) 
j j 

where in the second form we expressed [/ as a product of commuting unitary evo- 
lutions by moving the projectors into the exponential. A general unitary map can 
be thus be constructed from d propagators of the form exp{— iAj|(/)j)((/)j|}, one for 
each eigenvalue/eigenvector pair. These unitary propagators can now be constructed 
using state mappings. We begin by noting that there exists some Vj G SU (d) that 
satisfies 

^-iXj\cl^j)m ^ g-iA,-y;|0>{0|y,- ^ ytg-iA,|0>{0|y. 

J ^ ' 

where |0) is a fixed "fiducial state". The sole requirement on Vj is that |(0|l^|(/)j)p = 
1, i.e., it must be a mapping from to |0). Therefore, we can create the unitary 
propagator exp{—iXj\(/)j){(l)j\} by using a state preparation to map the eigenvector 
of [/, onto the fiducial state |0), applying the correct phase shift, and finally 
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mapping the fiducial state back to the eigenvector with the time-reversed state prepa- 
ration. A general unitary map is thus constructed via the sequence, 

U = yje-^^^l°><°ly^ . . . y2^e-^^^l°><°ly2n^e-^^^"^^"V^i. (5.3) 

Each of the propagators Vj is specified by a control waveform that generates a desired 
state mapping. One can efficiently find such control fields based on a numerical search 
that employs a simple gradient search algorithm, as described above. To generate 
an arbitrary element of SU{d)^ we require at most d such searches. Moreover, the 
full construction consists of 2d state preparations interleaved with d applications of 
the phase Hamiltonian, leading to an evolution that is only of order d times longer 
than a state mapping evolution. 

This construction places only two requirements on the Hamiltonian in addition to 
controllability. Firstly, the dynamics must be reversible such that if we can generate 
the unitary evolution Vj^ we can trivially generate the unitary Vj by time-reversing 
the control fields. Note that this is not the same as finding a state preparation that 
goes in the opposite direction, |0) there are many unitary propagators that 

map |0) so it is unlikely to find the unique operator Vj from a stochastic 

search. Generally, we can easily time reverse our controls if the Hamiltonian dynam- 
ics have no drift term. In some cases with a non-zero drift term it is still possible 
to time reverse controls, however, in this case we need to be able to find a rotating 
frame that removes the drift term while leaving the remaining Hamiltonian terms 
reversible. Secondly, we require access to a control Hamiltonian that applies an arbi- 
trary phase to one particular fiducial state |0) relative to all of the remaining states 
in the Hilbert space, exp{— iAj|0)(0|}. This latter requirement is the most restrictive, 
but can be implemented in a wide variety of systems. An example is discussed in 
5.2. 
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5.1.1 Subspace maps 

We have so far considered two kinds of maps on our rf-dimensional Hilbert space H: 
dxd unitary matrices and state-to-state maps. The former corresponds to a map U : 
TC ^ while the latter specifies a map on a one-dimensional space. Intermediate 
cases are also important. In particular, we are often interested in unitary maps that 
take subspace A of arbitrary dimension n to subspace according to T : A ^ B. 
Examples include the encoding of a logical qubit into a large dimensional Hilbert 
space {A ^ B) and a logical gate on encoded quantum information [A = B). Above 
we showed that the design of a fully-specified unitary matrix required search for d 
waveforms that define d state preparations (trivially a state mapping requires one 
such search). We show here how unitary maps on subspaces of dimension n can be 
constructed from exactly n such numerical solutions. 

Formally, a unitary map between two subspaces A and B of dimension n is defined 
as a map between between their orthonormal bases and {l&z)}, 

n 

Tn{A^B) = Yl ® (5.4) 

i=l 

where U± is an arbitrary map that preserves unitarity on the orthogonal complement 
Ai_ whose dimension is d—n. State preparation is the case n = 1; a full unitary matrix 
is specified when n — d. Clearly for n 7^ g? the map is not unique, with implications 
for the control landscape and the simplicity of numerical searches described above. 
As a first naive construction of T{A B), one might consider a sequence of one- 
dimensional state mappings, 

n 

T^{A^B) = \[T,{\ai)^\hi)). (5.5) 

This does not, however, yield the desired subspace map because each state mapping 
acts also on the orthogonal complement, so, e.g. \hi) is affected by Ti (|a2) I62)) 
and subsequent maps will move formerly correct basis vectors to arbitrary vectors 
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in the orthogonal component. We can resolve this problem by instead constructing 
subspace maps as a series all well-chosen rotations that maintain proper orthogonality 
conditions. 

To construct the necessary unitary operators, we make use of the tools described 
above: arbitrary state mapping based on an efficient waveform optimization and 
phase imprinting on a fiducial state. With these, we define the unitary map between 
unit vectors \a) and 

S (|a), \b)) = e-^"l^><^l = / - 2|(/))((/)|. (5.6) 

Here = N{\a) — where we have chosen the phases such that {b\a) is real and 
positive, and = 2(1 — {b\a)) is the normalization. This unitary operator has 

the following interpretation. In the two-dimensional subspace spanned by \a) and 
1 6), S' is a TT-rotation that maps S\a) = \b). In contrast to the state preparation map, 
Eq. (5.4) with n = 1, this map acts as the identity on the orthogonal complement to 
the space. This property is critical for the desired application. 

With these 2D primitives in hand, we can construct the subspace map according 
to the prescription, 

Tn{A ^ B) = Sn--- S2S1, (5.7) 

where Sk = S {\dk), \bk)) and 

\dj) = Sj_i . . . S2Si\aj). (5.8) 

This sequence does the job because each successive rotation leaves previously mapped 
basis vectors unchanged. To see this, we must show that at step j, the basis vec- 
tors I62), • • • , are unchanged by Sj. This will be true when this set is 
orthogonal to the vectors \aj) and \bj). Orthogonality to \bj) is trivial since the basis 
vectors of B are orthonormal. We must thus prove, {dj\bk) = 0, Vj > k. We can do 
this by induction. For an arbitrary fc, assume the conjecture is true for all j such 
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that jo > j > and thus Sj\bk) = \bk) up to j = jo. This imphes that {aj^^ilbk) = 
since, 



To complete our proof by induction, we must show that for any fc, the conjecture is 
true for j = k + 1. This foUows since, 

{o^k+i\bk) = {dk + 1|44 • • • ^k\bk) 



With this protocol we can construct unitary maps on a subspace of dimension n 
with optimized waveforms that corresponded to exactly n prescribed state prepara- 
tions. In the following section we apply these tools to qudit manipulations in atomic 
systems. 

5.2 Applications to the microwave rf system 

In this section, we apply our results to the control of the ground-electronic manifold of 
magnetic sublevels in alkali atoms discussed in chapter 3.2. In addition to an efficient 
method for designing and implementing state-to-state mappings, our protocol places 
certain requirements on the available control tools. Firstly, the system dynamics 
must be reversible so that we can trivially invert a state mapping. This is easily 
achieved through phase control. Secondly, we require phase imprinting on a single 
fiducial state. While this cannot be accomplished using solely microwave and rf- 
control, by introducing an excited electronic manifold, an off-resonant laser-induced 




(^jo+il4---44+i---4ol^^) 



(5.9) 



{ctk+i\cik) = 0. 



(5.10) 
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Figure 5.1: The hyperfine structure of ^^^Cs in the 6S1/2 ground state. Microwaves 
(blue) and rf magnetic fields (purple) provide controllable dynamics on the 16- 
dimensional Hilbert space. A detuned laser light shift (red) can be used to create 
a relative phase between the F = 4 and F = 3 manifolds. By considering controls 
on the subspace of the orange states we recover a system that satisfies the criteria 
proposed in chapter5.1. 



light-shift can achieve this goal. We restrict our system to one spin manifold (here 
the F = 3, but in principle either will do) and a single state from F = 4 manifold, e.g. 
|F = 4,m = 4), which acts as the fiducial state. By detuning far compared to the 
excited state line width of 5 MHz, but close compared to the ground-state hyperfine 
splitting of 10 GHz, we imprint a light shift solely on the |F = 4, m = 4) state with 
negligible decoherence. Using rf-magnetic fields to perform rotations in the F = 3 
manifold, and microwaves to couple to the fiducial state, we obtain controllable and 
reversible dynamics. Note that we may include the fiducial state in our Hilbert 
space, for a total of 8 sublevels, or treat it solely as an auxiliary state and restrict 
the Hilbert space to the 7-dimensional F = 3 manifold. 
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5.2.1 Constructing qudit unitary gates 

The standard paradigm for quantum information employs two-level systems - qubits 
- in order to implement binary quantum-logic based on SU (2) transformations. Ex- 
tensions beyond binary encodings in > 2 system - qudits - based of SU (d) trans- 
formations have also been studied and may yield advantages in some circumstances 
[68, 15, 69]. Of particular importance for fault-tolerant operation is implementation 
of these transformations through a finite set of "universal gates". Our goal here is 
to show how important members of the universal gate set can be implemented using 
our protocol. 

In choosing a universal gate set appropriate for error correction, it is natural to 
consider generalizations of the Pauli matrices X and Z which generate SU{2). The 
generalized discrete Pauli operators for SU (d) are defined 



Here © refers to addition modulo d and uu is the primitive dth root of unity, cu = 
exp{i27r/rf}. By considering the commutation relation of X and Z, the remaining 
generalized Pauli operators have the form cu^X^ Z^^ defining the elements of Pauli 
group for one qudit (up to a phase). This group is a discrete (finite dimensional) 
generalization of the Weyl-Heisenberg group of displacements on phase space. 

Another important group of unitary matrices in the theory of quantum error 
correction is the single qudit CliflFord group, given its relationship to stabilizer codes 
[68]. These group elements map the Pauli group back to itself under conjugation. 
Expressed in terms of their conjugacy action on X and Z, the generators of the 



m 

z\j) 



\j © 1) 



(5.11) 
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Clifford group for single qudits are 

HXH^ = Z, HZH^ = X-^ (5.12) 

SXS^ = XZ, SZS^ = Z (5.13) 
GaXGl = X^ G^ZGl = 

when gcd(a,(i) = 1 (5.14) 

H and S are direct generalization of the Hadamard and phase-gates famihar for 
qubits [ ]. The rf-dimensional H is the discrete Fourier transform 

and is a nonhnear phase gate 

S\j) = u;^U-m\j) ^-odd, (5.16) 
S\j)=u;^"%) J even. (5.17) 



The operator Ga is a scalar multiplication operator with no analog in the standard 
Clifford group on qubits, defined by 

Ga\j) = \aj), (5.18) 

where the multiplication is modulo d. The only such multiplication operator for 
2-level systems is the identity operator. 

While both the generalized Pauli and Clifford groups have utility in quantum com- 
puting, it is clear from their descriptions that unlike their qubit SU{2) counterparts, 
these unitary matrices do not arise naturally as the time evolution operators governed 
by typical Hamiltonians. This fact is not relevant to our unitary construction, which 
requires only knowledge of the operators' eigenvectors and eigenvalues. Using the 
time-dependent Hamiltonian dynamics with couplings illustrated in Fig. 5.1 we have 
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Figure 5.2: Optimized control fields for implementing the 7-dimensional Fourier 
transform on the F = 3 hyperfine manifold in ^^^Cs. The duration of the pulse 
is less than 1.2 ms and yields a unitary map that has an overlap of 0.9854 with the 
desired target. As an example, we show the action of the resulting unitary on the 
Z-eigenstates of angular momentum. The conjugate variable of is the azimuthal 
angle (j). If we Fourier transform a Z-eigenstate, a state with a well defined value of 
Fz^ we obtain a state that has a well defined value of 0, a squeezed state. 



engineered control fields to create the generators of both the Pauli and Clifford groups 
acting on the 7-dimension F = 3 hyperfine manifold. The duration of waveforms is 
approximately 1.5 ms, which is significantly shorter than the decoherence time of the 
system. In principle, the durations of these waveforms could be decreased by an order 
of magnitude or more by using more powerful control fields. Our objective function 
for creating a desired unitary W is the trace distance J\W] = Tr (W^U)^ where U 
is the unitary matrix generated by our control waveforms. Based on our protocol, 
employing state mappings that have fidelities of 0.99, our construction yields uni- 
tary maps that reach their targets with fidelities of J[Z] = 0.9866, J[X] = 0.9872, 
J[H] = 0.9854, J[S] = 0.9892 and J[Gs] = 0.9801. 

As an example, in Fig. 5.2 we show the control sequence for the discrete Fourier 
transform. The unitary map generated by this sequence should act to transform 
eigenstates of Z into eigenstates of X and vice versa. We illustrate this through a 
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Wigner function representation on sphere [ ]. The Z eigenstates are the standard 
basis of magnetic sublevels, whose Wigner functions are concentrated at discrete 
latitudes on the sphere, Fig. 5.2a. Applying the control fields to each of these states 
yields the conjugate states, with Wigner functions shown in Fig. 5.2b. These have 
the expected form. They are spin squeezed states concentrated at discrete longitudes 
conjugate to the Z eigenstates. The Z and X eigenstates are analogous to the number 
and phase eigenstates of the harmonic oscillator in infinite dimensions. 

5.2.2 Error-correcting a qubit embedded in a qudit 

The ability to generate unitary transformations on two-dimensional subspaces allows 
us to encode and manipulate a qubit in a higher dimensional Hilbert space in order 
to protect it from errors. Such protection can take a passive form through the 
choice of a decoherence-free subspace [71, 72], or active error correction through an 
encoding in a logical subspace chosen to allow for syndrome diagnosis and reversal 
[73, 8]. Typically, error protection schemes involve multiple subsystems (e.g. multiple 
physical qubits) to provide the logical subspace. While tensor product Hilbert spaces 
are generally necessary to correct for all errors under reasonable noise models, for a 
limited error model, one can protect a qubit by encoding it an a higher dimensional 
qudit [ ]. We consider such a protocol as an illustration of our subspace-mapping 
procedure. 

As an example, we consider encoding a qubit in the ground-electronic hyperfine 
manifold of ^^^Cs and protecting it from dephasing due to fiuctuations in external 
magnetic fields. In the presence of a strong bias in the z-direction, the spins are 
most sensitive to fiuctuations along that axis. For hyperfine qubits, one solution is 
to choose the bias such that two magnetic sublevels see no Zeeman shift to first order 
in the field strength (a "clock transition"). An alternative is to employ an active 
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error correction protocol analogous to the familiar phase-flip code [70]. 

We take our "physical qubit" computational basis to be the stretched states, |0) = 
|3,3^) and |1) = |4,4^), states easily prepared via optical pumping and controlled 
via microwave- drive rotations on the Bloch sphere Fig. 5.3(i). Here we have used the 
shorthand labeling the two quantum numbers |F, m^), and have denoted the relevant 
quantization axis by the subscript on the magnetic sublevel. Such states, however, 
are very sensitive to dephasing by fluctuations along the bias magnetic fleld, and such 
errors are not correctable. We choose as our encoded qubit basis stretched states 
along a quantization axis perpendicular to the bias (x-axis), {|0) = |3,3a^),|I) = 
|3, —3a;)}, Fig. 5.3(ii). Choosing this basis, a dephasing error in the z-direction acts 
to transfer probability amplitude into an orthogonal subspace. Such errors that can 
be detected and reversed without loss of coherence. 

Our error correction protocol works as follows. Consider an encoded qubit \ip) = 
a\0) + f3\l) . The error operator due to B-fleld fluctuations is the generator of rota- 
tions, Fz. Assuming a small rotation angle 26, when such an error occurs, the state 
of our encoded qubit is mapped to 

^\xP)+e («|3, 2,) + /?|3, -2,)) . (5.19) 

The error acts to spread our qubit between two orthogonal subspaces, \mx\ = 3 and 
\mx\ = 2, Fig. 5.3(iii). To diagnose the syndrome we must measure the subspace 
without measuring qubit. We can achieve this by coherently mapping the error 
subspace to the upper hyperflne manifold. Fig. 5.3(iv), followed by a measurement 
that distinguishes the two hyperflne manifolds, F = 3 and F = 4. Such a coherent 
mapping cannot be achieved through simple microwave-driven transitions since the 
bias fleld is along the z-direction while the encoded states are magnetic sublevels 
along the x-direction. We can instead use the construction of unitary operators on 
a subspace described in chapter 5.1 to design vr-rotations that take the error states 
to the upper manifold. This is tricky for our implementation because our protocol 
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Figure 5.3: (A) A schematic of the error 
correction protocol we have designed us- 
ing subspace maps. We track the basis el- 
ements of our encoded subspace, here |0) 
is red and |1) is blue, via their populations 
in the x and z bases. The different config- 
urations are explained in the text. In (B) 
we examine the performance of the error 
correction. On the x-axis we have the an- 
gle of rotation in the z-direction due to 
the magnetic field error. On the y-axis is 
the fidelity between the initial and post 
error states, as average over pure states 
drawn from the Harr measure. The blue 
line shows the fidelity of the error cor- 
rected states and the green the fidelity if 
the state had simply stayed in the sub- 



space |4,4^), 



|3,3, 



only included one magnetic sublevel in the F = 4 manifold so as to ensure proper 
phase imprinting. The solution is to switch the auxiliary state in the upper manifold 
between two different subspace maps. First, we consider the control system where 
14,4^) is our auxiliary state and perform a vr-rotation that maps |3,2a^) to |4,4^), 
leaving the rest of the space invariant. Then employ control on the system where 
|4, — 4^) is the auxiliary state and map |3, — 2^;) to |4, — 4^), with the identity on 



83 



Chapter 5. EfEciently Constructing Arbitrary Unitary Maps on Qudits 

the remaining space. A QND measurement of F collapses the state to the initially 
encoded state when the measurement result is F = 3, or to the state a|4,4a^) + 
/?|4, —4a;), if we find F = 4, Fig. 5.3(v). In the final step of the protocol, if an error 
occurred, we conditionally move the error subspace back to the encoded subspace, 
which can be achieved through reverse maps of the sort described above. 

We simulate here the coherent steps in the error correction protocol. These are 
implemented through our efficient search technique to construct subspace maps for 
the sequences 

{|4,4,),|3,3,)} ^ {|3,3,),|3,-3,)} 
{|3,2,),|3,-2,)} ^ {|4,4,),|4,-4,)} 
{|4,4,),|4,-4,)} ^ {|3,3,)|3,-3,)} 

Each of these maps are achieved through a sequence of SU (2) vr-rotations on a two- 
dimensional subspace that leave the orthogonal subspaces invariant. Starting from 
numerical searches for state preparation maps that have fidelity greater than 0.99, 
we obtain subspace maps with comparable fidelities. The performance of this error 
correction procedure is shown in Fig. 5.3B. We plot the fidelity between the initial 
state and the post-error-corrected state, averaged over random initial pure states of 
the physical qubit, versus the magnitude of the error as described by the rotation 
angle induced the stray magnetic field. Even with imperfect subspace transforma- 
tions the error correction protocol is significantly more robust than free evolution. 
Of course, like all quantum error correction protocols, we assume here that the time 
necessary for diagnosing the syndrome and correcting an error is sufficiently shorter 
than the dephasing time, so that the implementation of error correction does not 
increase the error probability. 

In practice, the most challenging step in the error correction protocol in this 
atomic physics example is measurement of the syndrome. This requires addressing of 
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individual atoms and measuring the F quantum number in a manner that preserves 
coherence between magnetic sublevels. In principle, this can be achieved through 
a QND dispersive coupling between an atom and cavity mode that induces an F- 
dependent phase shift on the field that could be detected [ ]. Alternatively, F- 
dependent fiuorescence from a given atom would allow this code to perform "error 
detection", without correction. 
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Summary and Outlook 



In this thesis I have studied the general principles of quantum control of finite dimen- 
sional quantum systems and their application to the control of alkali atomic spins. 
In Ch. 2 I discussed the general complexity of open-loop quantum control, focusing 
on the two control tasks: state preparation and unitary construction. In particu- 
lar, I provided a pedagogical review of a sequence of papers from the Rabitz group, 
[47, 48, 49, 50, 51, 52], on the topic of control landscape topology. The conclusion 
from these papers is that the problem of state preparation has a landscape that is 
very favorable with respect to local searches for optimal controls, while the landscape 
for the problem of unitary construction is more much complex, and makes finding 
optimal controls with the same types of searches unfeasible. This implies that we 
must utilize smarter algorithms to find optimal controls when construct full unitary 
maps. 

The platform with which we explored these control protocols was the alkali atomic 
spin system discussed in Ch. 3. I described two independent control systems devel- 
oped for these atomic spins. In the first, 3.1, the control was achieved through 
applied time-dependent magnetic fields that give rise to a Zeeman interaction, which 
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together with an "always-on" nonhnear hght shift, provided for fuU controUabihty. 
The Hilbert space controUed with this system was the 7-dimensional F = 3 ground- 
state hyperfine manifold in ^^^Cs. The second control system contained no optical 
fields and instead used oscillating magnetic fields at both rf and microwave frequen- 
cies to control the entire 16-dimensional electronic ground state 3.2. This system is 
favorable for the types of control we considered in later sections due to both the large 
number of tunable parameters, which led to many controllable configurations and the 
simple geometric description of the constituent Hamiltonians as representations of 
5u(2). Additionally, this systems dynamics are essentially coherent, which was not 
the case with the previous system where the laser light shift induces decoherence in 
the form of spontaneous photon scattering. 

In Ch. 4 I discussed open-loop state preparation. In this control task, the goal is 
to map a known fiducial state to an arbitrary target with unit fidelity. I developed 
an algorithm for finding good control waveforms in Ch. 4.1, which is based on simple 
gradient search techniques. We utilized this algorithm to construct state prepara- 
tions for both of the control systems in Ch. 3. By employing a nonlinear light shift 
in conjunction with time- varying magnetic fields, an experimental implementation 
of the state-preparation for waveforms of duration of about 0.5ms yielded the target 
with a fidelity on the order of 0.8 - 0.9. The difference between the fidelity in the 
optimization and in the experiment can be traced to known quantities, such as the 
precision in the density-matrix reconstruction protocol , inhomogeneities in the laser 
field or a rotation of the final state due to a mismatch between the state preparation 
and reconstruction waveforms. In addition, we looked at preparing squeezed states 
using our state preparation algorithm in comparison to the adiabatic technique pro- 
posed in [ ]. Our method for state preparation was about five times faster, which 
meant that there was less decoherence due to photon scattering. The imprecision in 
our application of these complicated control waveforms, however, produced squeezing 
that was slightly smaller than what was seen using the adiabatic technique. 
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With the microwave and rf control fields we were able to numerically construct 
state preparation protocols from our algorithm that were about an order of magnitude 
faster (50/xs - 150/xs) even though they acted on a space about twice as large (see 
Ch. 4.3). The main reason for the speedup is that we can increase the strength of 
the effective nonlinearity without fear of increased decoherence rates. We looked 
at the performance of a variety of scenarios, restricting some control parameters 
by, e.g., fixing the amplitudes of the fields or the number of resonant microwaves 
frequencies. Under certain conditions, restricted control yielded better performance. 
We suspect this is a numerical issue related to the complexity of searching a large 
dimensional control space. These unrestricted control system should contain higher 
fidelity waveforms, but for realistic parameters also contain more local optima. There 
is simply a larger space to sample from and the moderate length of our searches was 
fixed independently of the size of the control space. 

With regard to unitary construction, in Ch. 5 I described a protocol that uti- 
lizes stochastic searches to construct state preparations, as opposed to stochastically 
searching for full unitary maps. The computational resources for this algorithm scale 
only polynomially with the dimension of the system's Hilbert space and the duration 
of the control waveforms also scale polynomially with d. This hybrid search technique 
places only very mild restrictions on the types of Hamiltonians with which our pro- 
tocol is applicable as opposed to the constraints from other geometric constructions. 
The conditions for applicability are that the system dynamics are controllable, we 
can time-reverse our control fields, and that we can imprint an arbitrary phase on a 
single fiducial state |0). With this system we can construct not only unitary matrices 
on the full Hilbert space, but also maps on subspaces. With a subspace of dimension 
n < rf, the difficulty of the numerical search as well as the duration of the control 
scales like n/d with respect that of the full unitary construction. As an example we 
looked at constructing unitary maps in the microwave and rf control system 5.2. The 
most restrictive constraint on our system is the necessity of imprinting a phase on a 
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single state. To achieve this we considered control restricted to the F = 3 manifold 
with 1 4, 4) as an auxiliary state. We used a laser light shift to apply a phase to the 
F = 4 manifold relative to the F = 3, which in this restricted control space acts to 
imprint a phase on the state |4,4). With this control system we constructed gener- 
ators for the generalized Pauli and Clifford groups on the F = 3 manifold, e.g. the 
discrete Fourier transform. Starting from physically reasonable state preparations of 
0.99 we obtained unitary constructions with maps with a target fidelity of around 
0.98. We have not performed any detailed error analysis, but this scaling suggests 
that the fidelity of full unitary maps may go as the square of the state preparation 
fidelity. Additionally we looked at constructing subspace maps to enable correction 
of errors due to unknown z-magnetic fields for a qubit encoded in this larger spin. 

There are some obvious extensions to the control systems in this thesis that are 
already being pursued by myself and others. The microwave and rf control system 
has the nice property that the independent terms in the Hamiltonian look like the 
generators of SU{2) in overlapping subspaces described by the two hyperfine spin 
manifolds, F = 3 and F = 4, as well as the 2D subspace spanning these manifolds 
that is coupled by the microwave interaction. Much of the power of qubit control 
comes from our understanding of rotations on the 2-sphere. In work with Brian 
Mischuck, we have shown that it is possible to create general unitary operators from 
SU{1Q) in this control system out of sequences of SU{2) rotations. We are currently 
working to make these individual rotations robust to detuning and amplitude errors 
using robust composite pulse design techniques from NMR control [ , ] in order to 
create robust SU{1Q) evolutions. 

From a more general perspective, in Ch. 2.2, we discussed the relative complexity 
for stochastically searching for control waveforms that generate state preparations 
versus a full dimensional unitary map. It takes resources polynomial in d to find 
good state preparation waveforms but resources exponential in d to find unitary con- 
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struct ions. While it is possible to determine some of the properties of the landscape 
analytically, to get the scaling, numerics are required. It would be nice to understand 
the transition point between polynomial and exponential scaling in these problems. 
A state preparation map put constraints on one row of a unitary matrix as opposed to 
constructing unitary maps where the entire matrix is specified. Subspace mappings 
lie somewhere in between. A reasonable hypothesis would be that the dimension of 
the subspace must have some intermediate scaling with the Hilbert space dimension 
in order to require exponential resources. Confirming this conjecture would be useful 
since it would imply that stochastic search techniques would scale efficiently with the 
problem of 2D subspace mapping. 

Beyond the uses for 2D subspace for encoding a qubit into a qudit, the ability to 
search for 2D subspace maps directly would greatly improve the applicability of the 
unitary construction technique in Ch. 5. Finding natural Hamiltonians of the form 
|0)(0|, to imprint a phase on the fiducial state, is fairly difficult. In most cases one 
has to proceed by coupling to an ancilla subspace that is not part of the Hilbert space 
being controlled, and then producing the effective Hamiltonian |0)(0| by tracing out 
the ancilla. If the primitive for the phase gate was of the form |0) (0| — 1 1) (1 1 no ancilla 
subspace is required. In fact, in the examples we have discussed there is already such 
a term in our Hamiltonian in the form of the microwave coupling. We have been 
able to show that it is trivial to extend our unitary protocol to work with this phase 
primitive and so all that remains is to determine whether one can efficiently search 
for 2D subspace mappings. 

An important tool in developing these sorts of quantum control protocols that 
I haven't discussed in much detail in this manuscript is measurement. Prior work 
on this subject was carried out in the collaboration between the Deutsch and Jessen 
groups in the PhD works of Andrew Silberfarb [ ] and Greg Smith [ ]. An im- 
portant goal for the near future is to extend their work, applied in the context of 
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magnetic field and nonlinear light shift, to the microwave and rf control system. I 
have collaborated with Carlos Riofrio, who has recently begun to work on this exten- 
sion. In the density matrix reconstruction procedure in [^to], the observables must be 
driven through dynamics, in the Heisenberg picture, so as to span an informationally 
complete set of measurements. Without an informationally complete set there will 
always exist some density matrices that cannot be reconstructed with unit fidelity. 
In a real physical system, there will additionally be errors in the measurement record 
due to a sensitivity to external fields as well as Gaussian noise associated with fi- 
nite measurement statistics. In this case it is important not only to sample from an 
informationally complete set, but to sample in some unbiased manner. Optimizing 
the dynamics so as to drive the observables uniformly through an informationally 
complete set has proved exceptionally difficult, and is just barely possible in a seven 
dimensional Hilbert space with reasonable computational resources. This is a prob- 
lem when we would like to consider a 16-dimensional system. In this case we have 
been looking at dynamics that correspond to sampling from pseudorandom unitary 
matrix distributions, such as those that arise from quantum chaotic maps [ ], in the 
hopes that random unitary evolution will provide measurement records that are suffi- 
cient for reconstruction without requiring a huge computational overhead to optimize 
the dynamics. 

With the tools in this thesis, state preparation and unitary construction, as well 
as the ability to perform density matrix reconstruction, we reach a level of control 
that allows us to explore new and interesting physics which, as a physicist, is a 
primary goal of developing quantum control techniques. With the AC-Stark shift 
control system, in [ ], the authors were able to use the state preparation techniques 
developed in [ ] and the density matrix reconstructions methods from [43, 44] to 
explore quantum chaos in the quantum-kicked top. In a 7-dimensional Hilbert space, 
which is deep within the quantum regime, it is possible to see features of the classical 
phase space if one can prepare and measure atomic spins. 
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The AC-Stark shift system has truly been a workhorse for exploring quantum 
control techniques [42, 44, 41, 40], and in the future, I expect that the microwave and 
rf magnetic field system should be able to fill a similar role. Currently, this control 
system is being constructed in Poul Jessen's lab at the University of Arizona. With 
state preparation and unitary construction we can hope to see more explorations of 
quantum chaos in this system, as well as perhaps studies of into many-body physics. 
The techniques in this dissertation should also be applicable to single atoms, such 
as atoms trapped in an optical lattice which is a well-known paradigm for quantum 
computing. 
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Appendix A 

Generalized Wigner function 
representation 

In dealing with high dimensional spin systems, it is useful to be able generate graph- 
ical representations of the quantum states which give some geometric intuition. The 
spin coherent state Wigner function representation introduced by Agarwal [ ] pro- 
vides a generalization of the standard Wigner function based on harmonic oscillator 
coherent states used to describe infinite dimensional systems. Given a spin J, the 
spin coherent state Wigner function is essentially a multipole representation on the 
sphere defined as, 

= EETr[prf (A.i) 

k m 

where Yq^\9^ (j)) are the spherical harmonics, and Tq^\j) are the irreducible spherical 
tensors given by 

(^) = y^T^E^^' m + q\k,q-J,m)\J,m + q){JM- (A.2) 

m 

For a given spin, the indices describing non-trivial irreducible tensors run from < 
k < 2J and —k<q<k. These plots are useful visualization tools because they 



94 



Appendix A. Generalized Wigner function representation 

capture the effect of geometric rotations on the quantum state. Two quantum states 
that differ solely by a SU(2) rotation will generate Wigner functions that also differ 
from each other by the same physical rotation. 

We seek to generalize this to the case of a tensor product space of two spins (here 
electron and nuclear), equivalent to the direct sum of two irreducible representations 
of SU(2) in the hyperfine subspaces, F and F^ We achieve this by considering the 
expanded set of tensors defined by 



The range of the indices is now \F — F'\ < k < F -\- F' and —k < q < k. One 
can easily show that for two spin manifolds, the set of operators {Tg^^(F, F), 



Tf\F,F'),Tf\F',F),Tf\F\F')} comprises a complete orthonormal operator ba- 



sis for the tensor product space. We again can map these operators to the spherical 
harmonics, and for each state get four spherical Wigner functions: one each for the 
F and F^ manifolds, and two for the coherences between manifolds. We label them 
Wf,f^ and Wpf^F'- By the Hemiticity of the density operator, Wf,f' and 

Wf'^f contain redundant information and are complex, so one need only consider the 
real and imaginary part of Wf,f'^ yielding four real functions. 

We scale the radii of the spheres over which the Wigner function is plotted. 
For the functions that describe a given hyperfine manifold, we let the radius of the 
sphere equal the population in the subspace, Tr(PppPp). In order to set the radii 
of the spheres corresponding to the coherences between the manifolds, we look at 
the sum of the singular values of the off-block component of the density matrix, 
X^m^ m|p|F', m')p. This allows for nonequal dimensions of the two sub- 



T^^\F,F') 




m 



(A.3) 
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space. Additionally, we scale these "coherence spheres" by the ratio of the real versus 
imaginary parts of Wigner function. The primary purpose of doing this is to be able 
to distinguish between pure superpositions and incoherent mixtures between the two 
manifolds. 

To gain some intuition, we show examples of different states and different repre- 
sentations. Figure (A.l) shows bar charts of the absolute values of the density matrix 
elements for the six states: \il^)ai = |4,4) and \ip)aii = (|4,4) + |3, — 3))V^ are spin 
coherent states and their superposition; l'^)^^, and \ip)bii are superpositions of spin 
squeezed states in the two manifolds along different quadratures; \ip)ci^ and \ip)cii are 
coherent superpositions vs. incoherent mixtures of a "cat state" (|3, 3) + |3, — 3))/\/2 
in one manifold, and a Dicke state |4, 0) in the other. The corresponding Wigner 
functions are shown in Fig. (A. 2). From these plots we make the following observa- 
tions. When restricted to a subspace corresponding to a given hyperfine manifold, 
the Wigner functions on the diagonal have the familiar forms of SU (2) Wigner func- 
tions, with the radius of the sphere determining the total population in that subspace. 
The off-diagonal Wigner functions show the effect of the coherences, had the entire 
Hilbert space been determined by an irreducible representation. This is clearly seen 
in Fig. (A.2aii), where the coherences are of the familiar form for a superposition of 
"north" and "south" pole spin coherent states. The effect of geometric rotation is 
exhibited in \ip)bi and \ip)bn- The bar charts do not indicate any similarity between 
the states, while the Wigner functions are clearly related by a 90 degree rotation. 
Finally, the difference between coherent superpositions and incoherent mixtures of 
states in the two manifolds is clearly seen in Fig. (A. 2c). 
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(c.i) ^ (cii) 



Figure A.l: Representations of states with bar charts of the absolute values of the 
density matrix elements, (a.i) is a spin coherent state \ip)ai = |4,4) and (a.ii) is a 
superposition two oppositely oriented spin coherent states, one for each of the two 
manifolds, \ip)aii = ^jd^, 4) + |3, —3). In (b.i, b.ii) we show the effects of rotations 
on a superposition of spin squeezed states, each determined as the ground state of 
— Fy in the respective irreducible manifold. Finally, in (c.i) we have a coherent 
superposition of the state |4, 0) and a cat state ;^(|3,3) + |3, — 3)) and in (cii) we 
have an incoherent mixture of those two states. 
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(a.i) 




9 



9 O 



(a.ii) 



# 






(b.i) 



(b.ii) 







(ci) 



(c.ii) 



Figure A. 2: Representations of the six states shown in Fig. (A.I) by the generahzed 
spherical Wigner functions. Each state is represented by four spheres. The spheres 
on the diagonal are the standard SU(2) Wigner functions in the F = 4 (upper 
diagonal) and F = 3 (lower diagonal) irreducible subspaces. The radius of these 
spheres, ranging from zero to one, determines the total population in that subspace. 
The off-diagonal spheres represent the coherences between the two subspaces 
text). 



see 
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Controllability code 



In this appendix I present the Mathmatica Code I used to determine whether a 
Hamiltonian system is controllable. This consists of two functions. The first note- 
book ''Algsize" computes the iterated commutators of of the initial Hamiltonians 
and determines whether whey span su(rf). The second notebook "Makebasis" cre- 
ates a two canonical bases for different representations of $u{d) and saves them to a 
file. The first representation is as irreducible operators on a rf-dimensional Hilbert 
space, while the second is the so-called adjoint representation, or the linear operators 
corresponding to the commutator action of our operators. 

Algsize^) 

Set Directory ["/Users/smerkel / research / alg_gen_mathmatica" ] ; 

optovec[x_, basis_] := Module[{d, i , jj , temp}, 
d = Dimensions [basis ] [ [ 1 ] ] ; 
temp = Table[0 , {i , (d^2-l)}]; 
For[jj =1, jj <= (d^2 - 1) , jj++, 

temp[[jj]] = Re[Tr[x. basis [[All, All, jj]]]]; 
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1 temp 

] 

AlgSize [ hamils_ ] := Module[{ inittime , dims, dh , dalg , dinit , 
zeroish , 

fnamel , fname2 , canbasis , admats , algbasis , 
baseproject , target , mm, nn , 

test , comaction , numtest , stabil , ttime , 111 , ttt , 
test2 , algsize } , 
inittime = SessionTime [ ] ; 

dims = Dimensions [ hamils ] ; 
dh = dims [ [ 1 ] ] ; 
dalg = dh^2 - 1; 
dinit = dims [ [ 3 ] ] ; 
zeroish = 10. ^(-10) ; 

fnamel = "canbasis" <> ToString [dh ] ; 
fname2 = "admats" <> ToString [dh ] ; 
canbasis = Get[fnamel]; 
admats = Get[fname2]; 

algbasis = Table [0, {i, dalg}, {j, dalg}]; 
baseproject = IdentityMatrix [ dalg ] ; 

target = 1; 
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For [ram = 1, rara <= dinit , rao:tf+, 

test = optovec[ h ami Is [ [ All , All , mm] ], canbasis]; 
If [ target > 1 , 

For[nn = 1 , nn < target , nn++, 

test = test — ( algbasis [ [ All , nn]].test) 
algbasis [ [ All , nn ] ] ; 



If[test.test > zeroish , 

algbasis [[ All , target]] = test / Sqrt [ test . test ] 
baseproject = baseproject — Outer [Times , 

algbasis [[All, target]] , algbasis 
All , target ] ] ] ; 

target++; 



For [ram = 2, rara <= dalg , rara-h+, 

If [rara = target , Break [ ] ] ; 

If [target = (dalg + 1), Break []] ; 

coraaction = Table [0, {i, dalg}, {j, dalg}]; 

For [ 111 = 1, 111 <= dalg , 111++, 
coraaction = 
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comaction + algbasis [ [ 111 , ram]] adraats [ [ All , 
All, 111]] ; 

]; 

For[nn = 1 , nn <= (mm — 1) , nn++, 
If [target = (dalg + 1), Break []] ; 
test = comaction . algbasis [[ All , nn ] ] ; 



test2 = baseproject . test ; 

If [ test 2 . test2 > zeroish , 

For[ttt = 1, ttt < target, ttt++, 

test = test — ( algbasis [[ All , ttt]]. test) 

algbasis [ [ All , ttt ] ] ; 
If [test, test < zeroish , Break[]] 

]; 

If[test.test > zeroish, 

algbasis [[ All , target]] = test / Sqrt[test. 
test ] ; 

baseproject = baseproject — Outer [ 

Times, algbasis [[ All , target]], algbasis [[ All , 

target ] ] ] ; 
target++; 
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numtest = Tr[baseproject * Transpose [baseproject]]; 
stabil = If [target < dalg , 0, 

If [numtest > lO^(-lO) ,1,0] 

]; 

algsize = target — 1; 

ttime = SessionTime [ ] — inittime ; 
{algsize , stabil , ttime} 

]; 

(^saving cannonical basis and comm actions^) 
Set Directory ["/Users/smerkel/research / alg_gen_mathmatica" ] ; 
d = 2 

canbasis = Table[0. , {ii , d}, {jj , d}, {kk, d"2 — 1}]; 
canine = 1; 

For [ii = l, ii<d, ii ++, 

For[jj = 1, jj < ii + 1, jj++, 

canbasis[[jj , jj , canine]] = 1/Sqrt[ii"2 + ii]; 

]; 

canbasis[[ii + 1, ii + 1, canine]] = — i i /Sqrt [ i i " 2 + 

ii]; 

caninc++; 

]; 

For [ii = l, ii<d, ii ++, 

For [ jj = ii + 1, jj < (d + 1) , jj++, 
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canbasis[[ii, jj, canine]] =1/ Sqrt [ 2 ] ; 
canbasis[[jj, ii, canine]] =1/ Sqrt [ 2 ] ; 
eanine++; 



For [ i i = 1 , 
For [jj = 

canbasis 
canbasis 
caninc++ 



i < d, ii++, 

i + 1, jj < (d + 1) , jj++, 
[ii, jj , canine]] =-I/Sqrt[2]; 
[ jj , ii , canine ] ] = I/Sqrt [2] ; 



admats = Table[0. , {ii , d^2 - 1}, {jj , d^2 - 1}, {kk , d 
^2 - 1}]; 

For[kk = 1, kk < (d^2 - 2) , kk++, 

For[ jj = kk + 1, jj < (d^2 - 1) , jj++, 
For[ ii = jj + 1, ii < (d^2) , ii++, 

temp = Re[Tr[— I (canbasis [[ All , 
All , kk ] ] . canbasis [ [ 

All, All, jj]] - canbasis [[ All , All, jj]]. 
canbasis [ [ All , All , 

kk ]]). canbasis [[ All , All, ii]] ]]; 

admats [[ii , jj , kk ] ] = temp; 
admats [ [ j j , kk , i i ] ] = temp ; 



104 



Appendix B. Controllability code 



admats [ [ kk , i i , j j ] ] = temp ; 

admats [ [ i i , kk , j j ] ] — —temp ; 

admats [ [ kk , j j , i i ] ] = —temp ; 

admats [ [ j j , ii , kk ] ] = —temp ; 



canbasis = SparseArray [ canbasis ] ; 
admats = SparseArray [ admats ] ; 
canbasis » "canbasis" <> ToString[d]; 
admats » "admats" <> ToString[d]; 
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State preparation code 



In this appendix I present my Matlab code for creating state preparations via gradient 
search. This is specificaUy written for use in the microwave and rf magnetic field 
control system. The code is broken into a number of files which are presented here 

The first file is the script, "opt_fid" which defines all the system parameters. The 
variables are stored in the data structure opt_params, which has components: 
init_state: initial state of system (usually |4, 4)) 
tot_time: time of total pulse length 

samp_rate: the sampling rate for out Schrodinger integrator 

mw_type: what type of control fields we use (see make_hamils_fields) 

mw_amp: rabi frequency 

mw_slew: maximum microwave slew rate 

rf_amp: Larmor freq. 

rf_slew: maximum rf slew rate 

hamils: array of hamiltonians (see make_hamils_fields) 

varJnfo: this is something else that comes from make_hamils_fields. Basically three 
numbers [number of optimization variables, the number of variables that belong to 
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the rf fields, the number of resonant microwave frequencies], 
fields: control fields found from optimization 
target: target state, usually just use random states 
fid: fidelity of optimized field 

%opt.fid 

%this is my script to maximise the fidelity for target 
states 

%full phase and amplitude control for both microwaves and rf 

both X and 
%y coils for the rf . 

Wo%%%%%%%%%%%%%%%%%%% 

%just some useful folders and loading stuff 
load ( ' Units . mat ' ) ; 

%ahs.path = ^ / share / Seth/ mwrf. cluster ^ ;% only need to worry 

about this with 
% %cluster 
dat a_save_folder = ' fields / test ' ; 

m%%%%%%%%%%%%%%%%%%% 

%specifying target states and intitial state, 
Wo%%%%%%%%%%%%%%%%%%% 
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%each target state should he a collumn vector and should 

have associated 
%with it a name. 

%initial state, this is the state \4U> 
init_state = zeros(16,l); 
init_state(l,l) = 1; 

% Here are two ways to make the target state , either using 

random targets 
% or by using a previously created target library 

m%%%%%%%%%%%%%%%% 

% random targets 

m%%%%%%%%%%%%%%%% 

n_teststates = 2]% number of targret states 
target = zeros (16 , n_teststates ) ; 
targ_name = cell (1 , n_teststates ) ; 
for ii=l:n_teststates 

target (:,ii) = random_state ( 1 6) ; 

targ_name{ i i } = str cat (' rand int2str ( i i )) ; 

end ; 

Wo%%%%%%%%%%%%%%%%%%%%%%%%%% 
% predefined set of targets 

m%%%%%%%%%%%%%%%%%%%%%%%%%% 

%Can generate random library with make^state^lib .m. 
% load tar g .lib .mat ; 
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% 

% if size (target , 2) ^= length (targ^name) 
% error(^need to label targets 

% end; 



vo/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o^^^^ 
%sp ecifying optimization parameters 
Wo%%%%%%%%%%%%%%%%%%% 

%almost all the optimization parameters can he inputted as 

vectors . If you 
%put in vector the code will loop over those parameter 

values . 



%general parameters 

tot_time = [50] * us ; %^zme of total pulse length 

samp_rate = [0. 1] * us ; %the sampling rate for out Schrodinger 

integrator 
%not going to loop over samprate 
%samp-rate=[5]^us;%fake for test 

iters = 20]%For each set of parameters we run the 
optimizaion iters number usually around 20 

%of times and take the best value ( something else we don ^ t 
loop over) 



%nfiw parameters 

mw_type ={ '2rfap2struwap ' };%/ooA; up in make.hamils .m 
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mw_amp = [40] * kHz ; %ra6i freq for for stretched state 

transition 
mw_slew = [5]^ MS ]%microwave slew time 
%mw-slew = [25]^ us;%just to test 



%rf parameters 

rf_amp = [15]*kHz;%r/ Larmor frequency 
rf_slew = [10]*us;%r/ slew time 
%rf.slew = [25]^ us;%fake for test 



job_size = length ( tot _time ) * length (mw_type) * length (mw_amp) * 
length (mw_slew) . . . 
* length ( rf_amp ) * length (rf_slew)*size(target ,2)*iters; 

f inal_it er = 0; 

m%%%%%%%%%%%%%%%%%%% 

%now we loop over everything creating a data structure 
opt-params that 

%contains all the optimization paramters and saving it to a 
file 

Wo%%%%%%%%%%%%%%%%%%% 

for ttime = 1 : length ( tot _time ) 
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for mwt = 1 : length (mw_type) 

for mwa = 1 : length (mw_amp) 

for mws = 1 : length (mw.slew) 

for rfa = 1 : length ( rf_amp ) 

for rfs = 1 : length ( rf .slew ) 



m%%%%%%%%%%%%%%%%%%% 

%initi all zin g mag field parameters into temp file 

m%%%%%%%%%%%%%%%%%%% 

clear opt_params 

opt_params .init_state = init.state; 
opt_params .tot_time = tot_time(ttime); 
opt_params . samp_rate = samp_rate; 
opt_params . mw_type = mw_type{mwt } ; 
opt_params .mw_amp = mw_amp(mwa) ; 
opt _params . mw_slew = mw_slew (mws) ; 
opt_params . rf_amp = rf_amp ( r f a ) ; 
opt_params . r f _sle w = r f _sle w ( r f s ) ; 

% making matrices for hamitonians and the vector var^info . 
hamils just 

% contains the hamiltonians in an array, var^info is three 
numbers 

%[number of optimization variables ^ the number of those 
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V ariablesthat 

% belong to the rf fields ^ the number of resonant microwave 
frequencies ] . 

[ h ami Is ,var_info] = make_hamils_fields (opt_params ,0) ; 

opt _params . hamils = hamils ; 
opt_params . var _inf o = var_info ; 

%opt-params . fields is where we will store the optimized 

waveforms ^ for now 
%zero 

opt_params . fields = zeros(4+2*opt_params . var _info (3) 
l+opt_params . tot_time/ opt_params . samp_rate ) ; 

for targ = 1 : size ( target , 2 ) 

opt_params . target = target (: , targ) ; 
opt_params . f id = 0; 

final_iter = final_it er a counter to see how many 
fields we^ll be finding 

m%%%%%%%%%%%%%%%%%%% 

%making file names 

m%%%%%%%%%%%%%%%%%%% 
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%we give each set of parameters a unique name 

waveform.name = strcat ( t arg .name {tar g }, ,int2str(tot_time 

( ttime ) /us ) , . . . 
' us_ ' ,mw_type(mwt) , ' _ ' , int2str (mw_amp(mwa) /kHz) , 'kHz ' , 

int2str (10*mw_slew (mws) / us ) , . . . 
'(usdlO) ' ,int2str (rf_amp(rfa)/kHz) , 'kHz' , int2str ( r f _sle w ( r f s 

)/us) 
'us'); 

%final-f names tells us the filenames the fields will be 

stored in when the 
%optimization is finished . It will look to see if you 

already have a file 
%with that name so as not to erase stuff if you want to try 

the whole thing 
%a second time 

final_fnames{ final_iter } = strcat ( data_save_folder , V ' ^ 

waveform_name , ' . mat ' ) ; 
if exist(char(final_fnames{final_iter}),'file') ^= 2 
save (char(final_fnames{final_iter}) , ' opt_params ' ) ; 

end 



for ii =1: iters 
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m%%%%%%%%%%%%%%%%%%% 

%temp file names 

m%%%%%%%%%%%%%%%%%%% 

% these are temporary filenames tthat we use as an input to 
our optimizer . 

% The optimizer will delete these files after it ^s don with 
them . 

waveform_name = strcat ( targ_name{targ } , '_@ ' , int 2str ( t ot _t ime 
( ttime ) /us ) , . . . 
' us_ ' , mw_type (mwt) , ' - ' , int 2str (mw_amp(mwa) /kHz) , 'kHz ' , 

int2str ( 10* mw_slew (mws) / us) , . . . 
'(usdlO) ' ,int2str(rf_amp(rfa)/kHz) , 'kHz' , int 2str ( rf .slew 

( rfs )/us) , . . . 
'us ' ,int2str( ii )) ; 

all_fnames{final_iter ,ii} = strcat (data_save_folder,'/', 
waveform_name , ' . mat ' ) ; 

save (char(all_fnames{final_iter , ii}) , ' opt_params ' ) ; 
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end 
end 
end 
end 
end 
end 
end 
end 

m%%%%%%%%%%%%%%%%%%% 

%Optimizing 
Wo%%%%%%%%%%%%%%%%%%% 

% This is the part where all the optimization actually 

happens. There ^ s two 
% chunks of code. The un commented hit is for running it on 

a single machine 
% and the commented hit is how I was running things on our 

cluster. Since each 
% function call is independent it is a hig speed up to 

parralelize this code. 

all_fnames_no_del = all_fnames ; % not used for anything 
except my own diagnostics 
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Wo%%%%%%%%%%%%%%%%%%% 
% on a local machine 
Wo%%%%%%%%%%%%%%%%%%% 



for kk = 1 : it er s 

for ii = 1 : size ( all_fnames , 1 ) 

results = make_optim ( char ( all_fnames { ii ,kk}) ,char( 
final_fnames{ii })) ; 

end 
end 



% cluster code 

m%%%%%%%%%%%%%%%%%%% 
% 

% sched = findResource ( ^ scheduler \ Hype \ ^jobmanager ^ ) ; 
% %sehed — findResource ( ^ scheduler \ Hype \ Hocal ^);%more 
diagnostics 

% %can run paralel 

stuff locally 

% 
% 

% fid_max = zeros ( 1 ^ size ( all_f names , 1) ) ; 
% all. fnames.no. del = all.fnames ; 
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% 
% 
% 

% j = create Job ( sched ) ; 

% 

% p = {ahs-path , strcat ( ahs-path , ^ / \ d at a^save^f older ) }; 

% set(j, ^ PathDependencies \ p) ; 

% for kk = 1 : iters 

% for ii — 1 : size ( all-f names , 1 ) 

% createTask (j ^ @make.optim2 , 2, {char ( all. fnames{ii ^ kk} ) 

, char ( final-fnames{ H})}) ; 
% end 
% end 
% 
% 

% submit (j); 
% 

% waitForState ( j ) ; 

% results = get AllOutput Arguments ( j ) 

% 

% 

% destroy(j); 



The function "make_hamils_fields" contains basically all the physics of the prob- 
lem. If fflag 0, this function creates the Hamiltonians, as well as some variables 
describing the number of optimization variables. It does this based on the the physi- 
cal setup, which I label with "opt_params.mwtype" . If fflag ==1, this code will take 
some raw optimization variables and fit them with cubic splines to create physical 
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waveforms with the proper slew rates. This is also cased out by "mwtype". 

function [ hamils .fields ,var_info] = make_hamils_fields ( 

opt_params , fflag ,x) 
%this is the function in which I basically put all of the 

physics of the 

%prohlem . Basically , when I want to change the physical con( 
that is by 

%making a new opt.params . mwtype) I only have to change thing 
here . I use 

%this function for two different things depending on the 
fflcig . 

%If fflag is 0, this function makes the hamiltonians in an 
array called 

%hamils and stores some of the relevant variable about the 

optimization in 
%the vector var.info . 

%If fflag is 1, this makes control fields out of the 

optimization variable 
%x using cubic splines. 

if isstruct ( opt_params )==0 
load ( opt_params ) ; 

end 

mwtype = opt_params . mw_type ; 
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switch mwtype 

case ' 2rfap2struwap ' 

%two spatial rf directions with amplitude and phase 
control 

%microwaves ressonant on both stretched state transition 

amplitude and 
%phase control 

if fflag = 

all_mw_trans = [3 ,4; — 3 , — 4] ;%m_F for microwaves 
rel_amps = [1 ; 1] ; %s c aling factor for rabi 
frequencies . 

%Useful if you want freq other 
than stretched 

fup = 4; 
fdown = 3; 

grel = -1.00321; % just g4/g3 
ntrans = size ( all_mw_trans , 1 ) ; 

hamils -fields = zeros (2* ( fup+fdown + 1) ,2* ( fup+fdown 
+ 1) ,4+2*ntrans ) ; 

%rf hamiltonians 
upang = make_gen ( fup ) ; 
downang = make.gen ( fdown ) ; 
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hamils -fields (: ,1) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

hamils_fields (:,:,!) = [ upang . jx , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros (2* fdown + 1 ,2* fup + 1) , gr el * (downang . jx ) ] ; 
hamils .fields (: ,2) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

hamils -fields (: ,2) = [upang. jy , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros(2*fdown + l,2*fup + l),— grel *( downang. jy ) ] ; 
hamils -fields (: ,3) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

hamils -fields (: ,3) = [upang. jy , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros (2* fdown + 1, 2* fup + 1) , grel *( downang. jy ) ] ; 
hamils -fields (: ,4) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

hamils -fields (: ,4) = [upang. jx , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros (2* fdown + 1 ,2* fup + 1) grel * (downang . jx ) ] ; 



%nw hamiltonians 

for i i = 1 : ntrans ; 

[ mw JK , mw_y ] = uw_maker_int ( all_mw_trans ( ii ,:)); 
hamils -fields (:,: ,4 + 2* i i —1) = reLamps ( i i ) *mwjK 
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hamils -fields (:,: ,4 + 2* i i ) = rel_amps ( i i ) *mw_y ; 

end 

%something for the distribution of variable for this 
optimization 

%let ^s say var^info is a vector with components ( 

total number of variables 
%needed , number aloocated to rf fields, number of 

microwave transition ) 

nrf_vars = 4*(ceil (opt_params . tot_time/ opt_params . 

rf _sle w ) — 1) ; 
nmw_vars = 2* ntrans * ( ceil ( opt_params . t ot _t ime / 

opt_params . mw_slew) — 1); 
var_info = [ nr f _var s+nmw_vars , nrf _vars , ntrans ] ; 

elseif fflag==l 

ntrans=opt_params . var_info (3) ; 
nrf.vars = opt_params . var _info (2) ; 

hamils -fields = zeros (4+2* ntrans ,l + opt_params . 
tot _time / opt_params . samp_rate ) ; 
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rf_vars = reshape (x(l: nrf_vars) ,4, nrf.vars / 4) ; 
mw_vars = reshape (x (( nrf_vars +1) : end) ,2* ntrans , ( 
length (x) — nrf.vars) / (2*ntrans)) ; 

t = 0: opt.params . samp.rate : opt.params . tot.time ; 

r f t = : opt_params . tot_time/(l+nrf_vars / 4) : 
opt_params . tot_time ; 

mwt = 0: opt.params . tot.time/ (l + (length (x)— nrf _vars ) 



(2* ntrans ) ) : opt _params . tot _time ; 
for ii =1:2 

rf_mags = opt_params . rf_amp* rf _vars (2* ii — 1 ,:) ; 
rf_thets= 2*pi*cumsum( rf _vars (2* ii ,:)); 
rfin = rf_mags . * cos ( rf _t het s ) ; 
rfout = rf _mags . * sin ( rf _t het s ) ; 

hamils_fields(2*ii — 1,:) = spline(rft ,[0,rfin ,0] , 



hamils_fields(2*ii ,:) = spline(rft ,[0, rfout ,0] ,t 



end 

for i i = 1 : ntrans 

mw_mags = opt _params . mw_amp*mw_vars (2* ii — 1,:); 



/••• 



t); 
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mw_thets= 2*pi*ciimsum(mw_vars (2* ii ,:)); 
mwin = mw_mags . * cos ( mw_thets ) ; 
mwout = mw_mags . * sin ( mw_thets ) ; 

hamils -fields (4+2* ii — 1 , :) = spline (mwt, [0 , mwin 
,0] ,t) ; 

hamils _f ields (4+2* ii ,:) = spline (mwt , [0 , mwout 
,0] ,t); 



end 
end 

case '2rfa2struwa ' 

%two spatial rf directions with amplitude control 

and fixed phase 
%microwaves ressonant on both stretched state 

transition amplitude 
%control and fixed phase 
if fflag = 

all_mw_trans = [3,4;— 3,— 4]; 
rel_amps = [1 ]l]]%s caling factor for rabi 
frequencies 

fup = 4; 
fdown = 3; 

grel = -1.00321; % just g4/g3 
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ntrans = size ( all_mw_trans , 1 ) ; 

hamils -fields = zeros (2* ( fup+fdown + 1) ,2* ( fup+fdown 
+ 1) ,2+ntrans ) ; 

%rf hamiltonians 
upang = make_gen ( fup ) ; 
downang = make_gen ( fdown ) ; 

hamils -fields (: ,1) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

hamils-fields (:,:,!) = [ upang . jx , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros (2* fdown + 1, 2* fup + 1) , grel *( downang. jx) ] ; 
hamils -fields (: ,2) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

hamils -fields (: ,2) = [upang. jy , zeros (2* fup + 1 ,2* 

fdown + 1) ; 

zeros (2* fdown + 1 ,2* fup + 1) , gr el * (downang . jy ) ] ; 
%Tiw hamiltonians 

for i i = 1 : ntrans ; 

[ mw JK , mw-y ] uw-makei-int ( all-UiW-t rans ( ii ,:)); 
h ami Is -f ields ( : , : , 2 + i i —1) = rel-amps ( i i ) *mwjK ; 

end 
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%something for the distribution of variable for this 
optimization 

%let ^s say var.info is a vector with components ( 

total number of variables 
%needed , number aloocated to rf fields , number of 

microwave transition ) 

nrf_vars = 2*(ceil (opt_params . tot_time/ opt_params . 

rf _slew ) — 1) ; 
nmw.vars = ntrans *( ceil ( opt _params . tot _time / 

opt_params . mw_slew) — 1); 
var.info = [ nr f _var s+nmw_vars , nrf _vars , ntrans ] ; 



elseif fflag==l 

ntrans=opt_params . var_info (3) ; 
nrf_vars = opt_params . var _info (2) ; 

hamils -fields = zeros(2+ntrans ,l + opt_params . tot_time 
/ opt_params . samp_rate ) ; 

rf_vars = reshape (x(l: nrf_vars) ,2 , nrf_vars / 2) ; 
mw_vars = reshape (x (( nrf _vars +1) : end) , ntrans ,( length 
(x) — nrf_vars) / ntrans) ; 
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t = 0: opt_params . samp_rate : opt_params . tot_time ; 

rft = 0: opt_params . tot_time/(l+nrf_vars / 2) : 
opt_params . tot_time ; 

mwt = : opt_params . tot _time / (l + ( length (x) — nrf_vars) 

/••• 

( ntrans ) ) : opt_params . tot _time ; 

for ii =1:2 

rf_mags = opt _params . rf_amp* r f _var s ( ii ,:) ; 
hamils_fields(ii ,:) = spline(rft ,[0, rf_mags ,0] , t 

); 



end 



for i i = 1 : ntrans 

mw_mags opt_params . mw_amp*mw_vars ( ii ,:) ; 
hamils .fields (2+ ii — 1 , :) = spline (mwt , [0 , mw_mags 
,0] ,t); 



end 



end 

case '2rfp2struwp ' 
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%two spatial rf directions with fixed amplitude and 
phase control 

%microwaves ressonant on both stretched state transition 

fixed amplitude and 
%phase control 

if fflag = 

all_mw_trans = [3,4;— 3,— 4]; 
rel_amps = [1 \l]\%s caling factor for rabi 
frequencies 

fup 4; 
fdown = 3; 

grel = -1.00321; % just g4/g3 
ntrans = size ( all_mw_trans , 1 ) ; 

hamils .fields = zeros (2* ( fup+fdown + 1) ,2* ( fup+fdown 

+ 1) ,4+2*ntrans ) ; 

%rf hamiltonians 
upang = make.gen ( fup ) ; 
downang = make_gen ( fdown ) ; 

hamils -fields (: ,1) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

hamils_fields(:,:,l) = [ upang . jx , zeros (2* fup + 1 ,2* 
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fdown + 1) ; 

zeros (2* fdown + 1 ,2* fup + 1) , gr el * (downang . jx ) ] ; 
hamils -fields (: ,2) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

hamils -fields (: ,2) = [upang.jy , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros (2* fdown + 1 ,2* fup + 1) ,— grel * (downang . jy ) ] ; 
hamils -fields (: ,3) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

hamils -fields (: ,3) = [upang.jy , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros (2 *fdown + 1,2* fup + 1) , grel *( downang. jy ) ] ; 
hamils -fields (: ,4) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

liamils-fields(:,:,4) = [ upang . jx , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros (2* fdown + 1 ,2* fup + 1) grel * (downang . jx ) ] ; 

%nw hamiltonians 

for i i = 1 : ntrans ; 

[ mw JK , mw-y ] = uw-maker-int ( all-UiW-trans ( ii ,:)); 
hamils -fields (:,: ,4 + 2* i i —1) = reLamps ( i i ) *mwjK 

hamils -fields (:,: ,4 + 2* i i ) = rel-amps ( i i ) *mw-y ; 

end 
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%something for the distribution of variable for this 
optimization 

%let ^s say var^info is a vector with components ( 

total number of variables 
%needed , number aloocated to rf fields , number of 

microwave transition ) 

nrf_vars = 2*(ceil( opt_params . tot _time / opt_params . 

rf _slew ) — 1) ; 
nmw_vars = ntrans *( ceil ( opt _params . tot _time / 

opt_params . mw_slew) — 1); 
var_info = [ nr f _var s+nmw_vars , nrf _vars , ntrans ] ; 

elseif fflag==l 



ntrans=opt_params . var.info (3) ; 
nrf_vars = opt_params . var _info (2) ; 

hamils -fields = zeros(4+2*ntrans ,l + opt_params . 
tot _time / opt_params . samp_rate ) ; 

rf_vars = reshape (x(l: nrf_vars) ,2 , nrf_vars / 2) ; 
mw_vars = reshape (x (( nrf _vars +1) : end) , ntrans ,( length 
(x) — nrf.vars) /(ntrans)) ; 
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t = 0: opt_params . samp_rate : opt_params . tot_time ; 

rft = 0: opt_params . tot_time/(l+nrf_vars / 2) : 
opt_params . tot_time ; 

mwt = : opt_params . tot _time / (l + ( length (x) — nrf_vars) 

/••• 

( ntrans ) ) : opt_params . tot _time ; 

for ii =1:2 

rf_mags = opt _params . rf_amp ; 
rf_thets= 2>i<pi*ciimsum( rf _vars ( ii ,:)); 
rfin = rf.mags .* cos ( rf _thets ) ; 
rfout = rf_mags .* sin ( rf _thets ) ; 

hamils_fields(2*ii — 1,:) = spline(rft ,[0,rfin ,0] , 

t); 

hamils_fields(2*ii ,:) = spline(rft ,[0, rfout ,0] ,t 

); 

end 

for i i = 1 : ntrans 

mw_mags = opt _params . mw_amp ; 
mw_thets= 2*pi*cumsmn(mw_vars ( ii ,:)); 
mwin = mw_mags . * cos ( mw_thets ) ; 
mwout = mw_mags . * sin ( mw_thets ) ; 
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spline (mwt , [0 , mwout 



spline (mwt, [0 , mwin 



end 
end 

case ' 2rfaplstruwap ' 

%two spatial rf directions with amplitude and phase 



%microwaves ressonant on one stretched state transition 

amplitude and 
%phase control 

if fflag = 
all_mw_trans = [3,4]; 

rel_amps = [l]]%scaling factor for rabi frequencies 

fup = 4; 
fdown = 3; 

grel = -1.00321; % just g4/g3 
ntrans = size ( all_mw_trans , 1 ) ; 



control 
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hamils -fields = zeros (2* ( fup+fdown + 1) ,2* ( fup+fdown 
+ 1) ,4+2*ntrans ) ; 



%rf hamiltonians 
upang = make_gen ( fup ) ; 
downang = make_gen ( fdown ) ; 

hamils -fields (: ,1) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

hamils_fields (:,:,!) = [ upang . jx , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros (2* fdown + 1 ,2* fup + 1) , gr el * ( downang . jx ) ] ; 
hamils -fields (: ,2) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

hamils -fields (: ,2) = [upang. jy , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros (2* fdown + 1 ,2* fup + 1) gr el * (downang . jy ) ] ; 
hamils -fields (: ,3) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

hamils -fields (: ,3) = [upang. jy , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros (2* fdown + 1, 2* fup + 1) , grel *( downang. jy ) ] ; 
hamils -fields (: ,4) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

liamils-fields(:,:,4) = [ upang . jx , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros (2* fdown + 1 ,2* fup + 1) grel * (downang . jx ) ] ; 
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%fnw hamiltonians 

for i i = 1 : ntrans ; 

[ mw JK , mw_y ] = uw_maker_int ( all_mw_trans ( ii ,:)); 
hamils_fields (: ,4 + 2* ii —1) = rel_amps ( i i ) *mwjx: 

hamils -fields (:,: ,4 + 2* i i ) = rel_amps ( i i ) *mw_y ; 

end 

%something for the distribution of variable for this 
optimization 

%let ^s say var^info is a vector with components ( 

total number of variables 
%needed , number aloocated to rf fields, number of 

microwave transition ) 

nrf_vars = 4*(ceil (opt_params . tot_time/ opt_params . 

rf _slew ) — 1) ; 
nmw_vars = 2* ntrans *( ceil ( opt_params . tot _time / 

opt_params . mw_slew) — 1); 
var_info = [ nr f _var s+nmw_vars , nrf _vars , ntrans ] ; 

elseif fflag==l 
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ntrans=opt _params . var.info (3) ; 
nrf_vars = opt_params . var _info (2) ; 

hamils -fields = zeros(4+2*ntrans ,l + opt_params . 
tot _time / opt_params . samp_rate ) ; 

rf_vars = reshape (x(l: nrf_vars) ,4, nrf_vars / 4) ; 
mw_vars = reshape (x (( nrf_vars +1) : end) ,2* ntrans , ( 
length (x) — nrf.vars) / (2*ntrans)) ; 

t = 0: opt.params . samp.rate : opt.params . tot.time ; 

rft = 0: opt.params . tot_time/(l+nrf_vars / 4) : 
opt_params . tot_time ; 

mwt = 0: opt_params . tot_time/ (l + (length (x) — nrf _vars ) 

/••• 

(2*ntrans)) : opt_params . tot _time ; 
for ii =1:2 

rf_mags = opt_params . rf_amp* rf _vars (2* ii — 1 ,:) ; 
rf_thets= 2*pi*cunisum( rf _vars (2* ii ,:)); 
rfin = rf_mags . * cos ( rf _t het s ) ; 
rfout = rf_mags .* sin ( rf _thets ) ; 

hamils_fields(2*ii — 1,:) = spline(rft ,[0,rfin ,0] , 

t); 
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hamils_fields(2*ii ,:) = spline(rft ,[0,rfout ,0] ,t 

); 

end 

for i i = 1 : ntrans 

mw_mags opt _params . mw_amp*mw_vars (2* ii —1,:); 
mw_thets= 2*pi*cumsum(mw_vars (2* ii ,:)); 
mwin = mw_mags . * cos ( mw_thets ) ; 
mwout = mw_mags . * sin ( mw_thets ) ; 

hamils -fields (4+2* ii — 1 , :) = spline (mwt , [0 , mwin 
,0] ,t) ; 

hamils .fields (4+2* ii = spline (mwt, [0 , mwout 
,0] ,t); 



end 



end 



case '2rfalstruwa ' 

%two spatial rf directio\ 

and fixed phase 
%microwaves ressonant on 

transition amplitude 
%control and fixed phase 



with amplitude control 



one stretched state 
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if fflag = 
all_mw_trans = [3,4]; 

rel_amps = [l]]%scaling factor for rabi frequencies 

fup = 4; 
fdown = 3; 

grel = -1.00321; % just g4/g3 
ntrans = size ( all_mw_trans , 1 ) ; 

hamils -fields = zeros (2* ( fup+fdown + 1) ,2* ( fup+fdown 
+ 1) ,2+ntrans ) ; 

%rf hamiltonians 
upang = make_gen ( fup ) ; 
downang = make_gen ( fdown ) ; 

hamils -fields (: ,1) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

liamils_fields(:,:,l) = [ upang . jx , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros (2* fdown + 1, 2* fup + 1) , grel *( downang. jx) ] ; 
hamils -fields (: ,: ,2) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

hamils -fields (: ,: ,2) = [upang. jy , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros (2* fdown + 1 ,2* fup + 1) , grel *( downang. jy ) ] ; 
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%fnw hamiltonians 

for i i = 1 : ntrans ; 

[ mw JK , mw_y ] = uw_maker_int ( all_mw_trans ( ii ,:)); 
hamils -fields (:,:, 2 + i i —1) = rel_amps ( i i ) *mwjK ; 

end 

%something for the distribution of variable for this 
optimization 

%let ^s say var^info is a vector with components ( 

total number of variables 
%needed , number aloocated to rf fields, number of 

microwave transition ) 

nrf_vars = 2*(ceil( opt_params . tot _time / opt_params . 

rf _slew ) — 1) ; 
nmw_vars = ntrans *( ceil ( opt _params . tot _time / 

opt_params . mw_slew) — 1); 
var_info = [ nr f _var s+nmw_vars , nrf _vars , ntrans ] ; 

elseif fflag==l 

ntrans=opt_params . var_info (3) ; 
nrf_vars = opt_params . var _info (2) ; 
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hamils -fields = zeros (2+ntrans ,l + opt_params . tot_time 
/ opt_params . samp_rate ) ; 

rf_vars = reshape (x ( 1: nrf_vars) ,2,nrf_vars / 2) ; 
mw_vars = reshape (x (( nrf.vars +1) : end) , ntrans ,( length 
(x) — nrf_vars)/ntrans) ; 

t = 0: opt_params . samp_rate : opt_params . tot_time ; 

rft = 0: opt.params . tot_time/(l+nrf_vars / 2) : 
opt_params . tot_time ; 

mwt = : opt_params . tot _time / (l + ( length (x) — nrf_vars) 

/••• 

( ntrans ) ) : opt_params . tot _time ; 

for ii =1:2 

rf_mags = opt_params . rf_amp* r f _var s ( ii ; 
hamils_fields(ii ,:) = spline(rft ,[0, rf.mags ,0] , t 

); 

end 

for i i = 1 : ntrans 

mw_mags opt _params . mw_amp*mw_vars ( ii ,:) ; 
hamils _f ields (2+ ii — 1 , :) = spline (mwt , [0 , mw_mags 
,0] ,t); 
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end 
end 

case '2rfplstruwp ' 

%two spatial rf directions with fixed amplitude and 
phase control 

%microwaves ressonant on one stretched state transition 

fixed amplitude and 
%phase control 

if fflag = 
all_mw_trans = [3,4]; 

rel_amps = [l]]%scaling factor for rabi frequencies 

fup = 4; 
fdown = 3; 

grel = -1.00321; % just g4/g3 
ntrans = size ( all_mw_trans , 1 ) ; 

hamils -fields = zeros (2* ( fup+fdown + 1) ,2* ( fup+fdown 
+ 1) ,4+2*ntrans ) ; 

%rf hamiltonians 
upang = make_gen ( fup ) ; 
downang = make.gen ( fdown ) ; 
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hamils -fields (: ,1) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

hamils_fields (:,:,!) = [ upang . jx , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros (2* fdown + 1 ,2* fup + 1) , gr el * (downang . jx ) ] ; 
hamils .fields (: ,2) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

hamils -fields (: ,2) = [upang. jy , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros(2*fdown + l,2*fup + l),— grel *( downang. jy ) ] ; 
hamils -fields (: ,3) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

hamils -fields (: ,3) = [upang. jy , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros (2* fdown + 1, 2* fup + 1) , grel *( downang. jy ) ] ; 
hamils -fields (: ,4) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

hamils -fields (: ,4) = [upang. jx , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros (2* fdown + 1 ,2* fup + 1) grel * (downang . jx ) ] ; 



%nw hamiltonians 

for i i = 1 : ntrans ; 

[ mw JK , mw_y ] = uw_maker_int ( all_mw_trans ( ii ,:)); 
hamils -fields (:,: ,4 + 2* i i —1) = reLamps ( i i ) *mwjK 
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hamils -fields (:,: ,4 + 2* i i ) = rel_amps ( i i ) *mw_y ; 

end 

%something for the distribution of variable for this 
optimization 

%let ^s say var^info is a vector with components ( 

total number of variables 
%needed , number aloocated to rf fields, number of 

microwave transition ) 

nrf_vars = 2*(ceil (opt_params . tot_time/ opt_params . 

rf _sle w ) — 1) ; 
nmw_vars = ntrans *( ceil ( opt _params . tot _time / 

opt_params . mw_slew) — 1); 
var_info = [ nr f _var s+nmw_vars , nrf _vars , ntrans ] ; 

elseif fflag ==1 

ntrans=opt_params . var_info (3) ; 
nrf.vars = opt_params . var _info (2) ; 

hamils -fields = zeros (4+2* ntrans ,l + opt_params . 
tot _time / opt_params . samp_rate ) ; 
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rf_vars = reshape (x(l: nrf_vars) ,2 , nrf.vars / 2) ; 
mw_vars = reshape (x (( nrf_vars +1) : end) , ntrans ,( length 
(x) — nrf.vars) /(ntrans)) ; 

t = 0: opt.params . samp.rate : opt.params . tot.time ; 

r f t = : opt_params . tot_time/(l+nrf_vars / 2) : 
opt_params . tot_time ; 

mwt = 0: opt.params . tot.time/ (l + (length (x)— nrf _vars ) 



(ntrans ) ) : opt _params . tot _time ; 

for ii =1:2 

rf_mags = opt _params . rf_amp ; 
rf_thets= 2*pi*cumsum( rf _vars ( ii ,:)); 
rfin = rf_mags . * cos ( rf _t het s ) ; 
rfout = rf _mags . * sin ( r f _t het s ) ; 

hamils_fields(2*ii — 1,:) = spline(rft ,[0,rfin ,0] , 



hamils_fields(2*ii ,:) = spline(rft ,[0, rfout ,0] ,t 



/••• 



t); 



end 



for 



i i = 1 : ntrans 



mw_mags = opt _params . mw_amp ; 
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mw_thets= 2*pi*ciimsum(mw_vars ( ii ,:)); 
mwin = mw_mags . * cos ( mw_thets ) ; 
mwout = mw_mags . * sin ( mw_thets ) ; 




spline (mwt , [0 , mwout 



spline (mwt, [0 , mwin 



end 
end 

case '2rfaplstruw0 ' 

%two spatial rf directions with amplitude and phase 



%microwaves ressonant on both stretched state transition 

fixed amplitude and 
%fixed phase always on^^ 

if fflag = 

all_mw_trans = [3,4]; 

rel_amps = [l]]%scaling factor for rabi frequencies 

fup = 4; 
fdown = 3; 



control 
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grel = -1.00321; % just g4/g3 
ntrans = size ( all_mw_trans , 1 ) ; 

hamils -fields = zeros (2* ( fup+fdown + 1) ,2* ( fup+fdown 
+ 1) ,2+ntrans ) ; 

%rf hamiltonians 
upang = make_gen ( fup ) ; 
downang = make.gen ( fdown ) ; 

hamils -fields (: ,1) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

hamils_fields(:,:,l) = [ upang . jx , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros (2* fdown + 1, 2* fup + 1) , grel *( downang. jx) ] ; 
hamils -fields (: ,2) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

hamils -fields (: ,2) = [upang. jy , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros (2* fdown + 1 ,2* fup + 1) grel * (downang . jy ) ] ; 
hamils -fields (: ,3) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

hamils -fields (: ,3) = [upang. jy , zeros (2* fup + 1 ,2* 

fdown + 1) ; 

zeros (2* fdown + 1 ,2* fup + 1) , gr el * (downang . jy ) ] ; 
hamils -fields (: ,4) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 
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hamils -fields (: ,4) = [upang.jx , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros (2* fdown + 1 ,2* fup + 1) grel * (downang . jx ) ] ; 
%nw hamiltonians 

hamils -fields (:,:, 5 ) = uw_maker_int ( all_mw_trans 

(1,0); 

%something for the distribution of variable for this 

optimization 

%let ^s say var.info is a vector with components ( 

total number of variables 
%needed , number aloocated to rf fields , number of 

microwave transition ) 

nrf.vars = 4*(ceil (opt_params . tot.time/ opt_params . 

rf _sle w ) — 1) ; 
nmw.vars = 0; 

var_info = [ nr f _var s+nmw_vars , nrf _vars , ntrans ] ; 

elseif fflag ==1 

ntrans=opt_params . var_info (3) ; 
nrf_vars = opt_params . var _info (2) ; 



145 



Appendix C. State preparation code 



hamils -fields = zeros (5 , 1 + opt _params . t ot _t ime / 
opt_params . samp.rate) ; 

rf.vars = reshape (x ( 1: nrf.vars) ,4, nrf.vars / 4) ; 

t = 0: opt_params . samp_rate : opt_params . tot_time ; 

rft = 0: opt_params . tot_time/(l+nrf_vars / 4) : 
opt_params . tot_time ; 

for ii =1:2 

rf_mags = opt_params . rf_amp* rf_vars(2*ii — 1,:); 
rf_thets= 2>Kpi*cumsum( r f _vars (2* ii ,:)); 
rfin = rf_mags . * cos ( rf _t het s ) ; 
rfout = rf _mags . * sin ( r f _t het s ) ; 

hamils_fields(2*ii —1,:) = spline(rft ,[0,rfin ,0] , 



hamils_fields(2>Kii = spline(rft ,[0, rfout ,0] ,t 



hamils_fields (5 , : ) = opt_params .mw_amp*ones (1 ,1 + 



t); 



end 
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opt_params .tot_time/ opt_params . samp .rate ) ; 

end 

case '2rfalstruw0 ' 

%two spatial rf directions with amplitude control 
and fixed phase 
%microwaves ressonant on both stretched state transition 

fixed amplitude and 
%fixed phase ^^always on^^ 

if fflag = 

all_mw_trans = [3,4]; 
rel_amps = [l]]%scaling factor for rabi frequencies 

fup = 4; 
fdown = 3; 

grel = -1.00321; % just g4/g3 
ntrans = size ( all_mw_trans , 1 ) ; 

hamils -fields = zeros (2* ( fup+fdown + 1) ,2* ( fup+fdown 
+ 1) ,2+ntrans ) ; 

%rf hamiltonians 
upang = make_gen ( fup ) ; 
downang = make_gen ( fdown ) ; 
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hamils -fields (: ,1) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

hamils_fields (:,:,!) = [ upang . jx , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros (2* fdown + 1, 2* fup + 1) , grel *(downang. jx) ] ; 
hamils -fields (: ,2) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

hamils -fields (: ,2) = [upang. jy , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros (2* fdown + 1 ,2* fup + 1) , grel *(downang. jy ) ] ; 
%Tiw hamiltonians 

hamils -fields (:,:, 5 ) = uw_maker_int ( all_mw_tr ans 

(1,0); 

%something for the distribution of variable for this 
optimization 

%let ^s say var.info is a vector with components ( 

total number of variables 
%needed , number aloocated to rf fields^ number of 

microwave transition ) 

nrf_vars = 2*(ceil( opt_params . tot _time / opt_params . 

rf _sle w ) — 1) ; 
nmw_vars = 0; 

var_info = [ nr f _var s+nmw_vars , nrf _vars , ntrans ] ; 
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elseif fflag==l 

ntrans=opt_params . var.info (3) ; 
nrf_vars = opt_params . var _info (2) ; 

hamils .fields = zeros (5 , 1 + opt _params . t ot _t ime / 
opt_params . samp_rate) ; 

rf_vars = reshape (x (1: nrf_vars) ,2 , nrf.vars / 2) ; 

t = 0: opt.params . samp.rate : opt.params . tot.time ; 

rft = 0: opt.params . tot_time/(l+nrf_vars / 2) : 
opt_params . tot_time ; 

for ii =1:2 

rf_mags = opt_params . rf_amp* r f _var s ( ii ,:) ; 
hamils_fields(ii ,:) = spline(rft ,[0, rf_mags ,0] , t 

); 

end 

hamils_fields(5,:) = opt_params . mw_amp* ones (1,1 + 
opt_params .tot_time/ opt_params . samp _r ate ) ; 
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end 

case '2rfplstruw0 ' 

%two spatial rf directions with fixed amplitude and 
phase control 

%microwaves ressonant on both stretched state transition 

fixed amplitude and 
%fixed phase always on^^ 

if fflag = 

all_mw_trans = [3,4]; 
rel_amps = [l]]%scaling factor for rabi frequencies 

fup = 4; 
fdown = 3; 

grel = -1.00321; % just g4/g3 
ntrans = size ( all_mw_trans , 1 ) ; 

hamils -fields = zeros (2* ( fup+fdown + 1) ,2* ( fup+fdown 
+ 1) ,2+ntrans ) ; 

%rf hamiltonians 
upang = make_gen ( fup ) ; 
downang = make_gen ( fdown ) ; 

hamils .fields (: ,1) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 
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hamils_fields (:,:,!) = [ upang . jx , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros (2* fdown + 1 ,2* fup + 1) , grel *(downang. jx) ] ; 
hamils -fields (: ,2) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

hamils -fields (: ,2) = [upang. jy , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros (2* fdown + 1 ,2* fup + 1) grel * (downang . jy ) ] ; 
hamils -fields (: ,3) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

hamils -fields (: ,3) = [upang. jy , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros (2* fdown + 1 ,2* fup + 1) , gr el * ( downang . j y ) ] ; 
hamils -fields (: ,4) = zeros (2* ( fup+fdown + 1) ,2*(fup+ 
fdown + 1) ) ; 

liamils-fields(:,:,4) = [ upang . jx , zeros (2* fup + 1 ,2* 
fdown + 1) ; 

zeros (2* fdown + 1 ,2* fup + 1) grel * (downang . jx ) ] ; 
%7iw hamiltonians 

hamils -fields (:,:, 5 ) = uw_maker_int ( all_mw_trans 

(1,0); 

%something for the distribution of variable for this 
optimization 

%let ^s say var.info is a vector with components ( 
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total number of variables 
%needed , number aloocated to rf fields^ number of 
microwave transition ) 

nrf.vars = 2*(ceil( opt_params . tot _time / opt_params . 

rf _slew ) — 1) ; 
nmw_vars = 0; 

var_info = [ nr f _var s+nmw_vars , nrf _vars , ntrans ] ; 

elseif fflag==l 

ntrans=opt_params . var_info (3) ; 
nrf.vars = opt_params . var _info (2) ; 

hamils -fields = zeros (5 , 1 + opt _params . t ot _t ime / 
opt_params . samp_rate) ; 

rf_vars = reshape (x(l: nrf_vars) ,2 , nrf_vars / 2) ; 

t = 0: opt_params . samp_rate : opt_params . tot_time ; 

rft = 0: opt_params . tot_time/(l+nrf_vars / 2) : 
opt_params . tot_time ; 

for ii =1:2 
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rf_mags = opt _params . rf_amp ; 
rf_thets= 2*pi*cumsimi( rf _vars ( ii ,:)); 
rfin = rf.mags . * cos ( rf _t het s ) ; 
rfout = rf_mags .* sin ( rf _thets ) ; 

hamils_fields(2*ii — 1,:) = spline(rft ,[0,rfiii ,0] , 

t); 

hamils_fields(2*ii = spline(rft ,[0, rfout ,0] ,t 

); 

end 

hamils_fields(5,:) = opt_params . mw_amp* ones (1,1 + 
opt_params .tot_time/ opt_params . samp .rate ) ; 



end 
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end 
end 



function [ mw jk , mw_y ] = uw_maker_int (mwtran) 

%little function to make the pauli operators between the 

correct m^F states 
fup = 4; 
fdown = 3; 

mwjx = zeros (2* ( fup+fdown + 1) ,2* ( fup+fdown + 1) ) ; 
mwjK(fup + 1 + mwtran(2) , 2*fup + 1 + fdown + 1 + mwtran(l)) 
= 1/2; 

mwjx:(2*fup + 1 + fdown + l+mwtran(l) , fup + 1+mwtran ( 2 ) ) = 

1/2; 

mw_y = zeros (2* ( fup+fdown + 1) ,2* ( fup+fdown + 1) ) ; 

mw_y(fup + l+mwtran(2) , 2>Kfup + 1 + fdown + 1+mwtran ( 1 ) ) = — 

i/2; 

mw_y(2*fup + 1 + fdown + l+mwtran(l) , fup + 1+mwtran (2) ) = i 

/2; 



end 
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"make_optim" basically takes an input file, optimizes the control waveform with 
"fmincon", deletes the input file, and conditional on the new waveform being better 
than previous waveforms save it to the specified save file location. 

function [t_timing, best_fid] = make_optim ( fname , save_name ) 
%make optim basically does all the optimization. It takes 
an input file 

%from fname, finds an optimal state preparation and staores 
it to save. name 

%conditional one the new fidelity being higher than whatever 
was previously 

%in save^name. Function outputs the fidelity as well as the 

time it took 
%the program to run. 

init_time = cputime ; 

load ( save_name ) ; 

past.fid = opt_params . f id ; 

if past_fid > 0.99 




%I decided we ^ d never need a waveform with a fidelity 

higher than 0.99 
%so if we already have a good waveform from a previous 

optimization 
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%this function doesn^t run an optimization You can 

change this value to whatever you 
%want , hut if you We running a big batch of 

optimizations some will 
%find good fields before the others so you^d like to not 
waste 



%resources optimizing something that ^ s already pretty 
good . 




best_fid = past_fid; 

delete (fname) ;%remo^'e5 temp file 

t_timing = 0; 

else 

load (fname) ; 

fidforpsi = @(x) — fid_mwrf ( opt_params , x) ]%objective for 

optimization 
lb = — ones(l, opt_params . var.info (1) ) ; 
ub = ones (1 , opt_params . var_info (1) ) ; 

rand_vars = rand ( 1 , opt.params . var_info (1) ) ; %random seed 
vars_lmax=fmincon (fidforpsi , rand_vars , [ ] , [ ] , [ ] , [ ] , lb , ub 

,[] 

optimset ( 'TolX ' ,le-3, 'TolFun ' ,le-3, 'Display ' , ' iter ' ) ) 
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%varsAmax = r and. vars ;% just for diagnostics optimization 

takes a long time 
fields = make_hamils_fields (opt_params ,1 , vars_lmax) ]%make 

fields from optimum 
opt_params . fields = fields; 

opt_params .fid = fid_mwrf ( opt _params ) \ %calculate fidelity 
best_fid = opt_params . f id ; 

% is this better than the previous optimium? 
if best_fid > past_fid 

save (save_name , 'opt_params ') ; 

else 

best_fid = past_fid; 

end 

delete (fname) ; %remoi;e temp file 
t_timing = (cputime— init _time ) /60; 
end 

''fid_mwrf ' calculates the fidelity of a state preparation. This can be called either 
with the data structure opt_params or a file name as an input. 

function fid = fid_mwrf ( opt_params , x) 

%will output fidelity of state preparation. opt_params is 
the data 

%structure opt^params in the optimization , hut can also he 
the filename 

%where opt^params is stored. I call this with the filename 
input after I 
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%have optimized fields to see that everything checks out. 
During the 

%optimization , espcially on a cluster architecture ^ it is 

fairly exspensive 
%to load things over and over again. x is the optimization 

variables , and 

%can he left out if you We using this function outsid eof 
the optimization . 

if isstruct ( opt_params )==0 
load ( opt_params ) ; 

end 

if nargin > 1 

%makes fields out of the optimization variables and puts 

them in opt_params 
fields = make_liamils_fields (opt_params , 1 ,x) ; 
opt_params . fields = fields; 
end 

%calls the schrodinger evolution 
psi_f = unit _evol_mwrf ( opt_params ) ; 

%calculate fidelity 

fid = abs ( opt_params . target '*psi_f); 

"unit_evoLmwrf ' is a Schrodinger integrator, 
function psi_f = unit _evol_mwrf ( opt_params ) 
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%just a simple schrodinger integrator 

if isstruct ( opt_params )==0 
load ( opt_params ) ; 

end 

psi_f = opt_params . init _st ate ; 

for ii = 1 : size ( opt.params . fields ,2) 
ht =0; 

for jj =1: size ( opt_params . fields ,1) 

ht = ht + opt_params . fields (jj , i i ) . * opt_params . 
hamils (: , jj ) ; 

end 

psi_f = expm(— i *opt_params . samp_rate *ht ) * p s i _f ; 

end ; 

"make_gen" provides generator's of angular momentum on an arbitrary spin, 
function Anggen = make_gen(s) 

%generates the irrdeuciable angular momentum operators for a 

spin — s system 

d = 2*s + l; 
Anggen . jx=zeros (d) ; 
for m=l:d 
for n = l:d 

i f (mH-l==n ) 
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Anggen . jx (m, n) =(1/2) *sqrt ( (d^n) *m) ; 
Anggen . jx (n ,m) =(1/2) *sqrt ( (d^n) *m) ; 

end ; 
end ; 
end ; 

Anggen . jy=zeros (d) ; 
for m=l:d 
for n = l:d 

if (n^l=n) 

Anggen . j y (m, n)=— i * ( 1 / 2) * sqrt ( ( d^n) *m) ; 
Anggen . jy (n ,m)=i * (1/2) *sqrt ( (d^n) *m) ; 

end ; 
end ; 
end ; 

Anggen . j z=zeros (d) ; 
for m =0:(d-l) 

Anggen. jz (m+l,m+l) = (d-l)/2 - m; 

end ; 

clear m n d 
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